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The method of matched asymptotic expansions is used to study the canonical problem of 
steady laminar flow through a narrow two-dimensional channel blocked by a tight-fitting 
finite-length highly permeable porous obstacle. We investigate the behaviour of the local 
flow close to the interface between the single-phase and porous regions (governed by the 
incompressible Navier-Stokes and Darcy flow equations, respectively). We solve for the 
flow in these inner regions in the limits of low and high Reynolds number, facilitating an 
understanding of the nature of the transition from Poiseuille to plug to Poiseuille flow in 
each of these limits. Significant analytical progress is made in the high-Reynolds-number 
limit, and we explore in detail the rich boundary layer structure that occurs. We derive 
general results for the interfacial stress and for the conditions that couple the flow in the 
outer regions away from the interface. We consider the three-dimensional generalization 
to unsteady laminar flow through and around a tight-fitting highly permeable cylindrical 
porous obstacle within a Hele-Shaw cell. For the high-Reynolds-number limit, we give the 
coupling conditions and interfacial stress in terms of the outer flow variables, allowing 
information from a nonlinear three-dimensional problem to be obtained by solving a linear 
two-dimensional problem. Finally, we illustrate the utility of our analysis by considering 
the specific example of time-dependent forced far-held flow in a Hele-Shaw cell containing 
a porous cylinder with a circular cross-section. We determine the internal stress within 
the porous obstacle, which is key for tissue engineering applications, and the interfacial 
stress on the boundary of the porous obstacle, which has applications to biohlrn erosion. 
In the high-Reynolds-number limit, we demonstrate that the fluid inertia can result in 
the cylinder experiencing a time-independent net force, even when the far-held forcing is 
periodic with zero mean. 
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1. Introduction 

Flow through a narrow channel containing a porous blockage is a canonical problem 
with numerous practical applications. For example, this how configuration is used : to pu¬ 


rify colloids in separation science, where it is also k nown as dead-end filtration (jVan der Bruggen et al 


2003 : McCarthy et al. 1998 : Bessiere et al. 20051) ; to deliver nutrient to cells growing in 


a porous tissue co nstruct in tissue engineering ( El Hai et al. Il990t lO’Dea et c/J.1 120091 


Jaasma et all 120081) : and to erode porous biohlms that have grown within a pipe, where 
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the e r osion is dependent on t he interfacial shear stress ( Picioreanu et oiJl2001 : Duddu et a,l 
20091 : Telemann et al. 2004 1. In many of these problems, the flow near the interface be¬ 
tween the fluid and the porous blockage is of experimental or industrial interest. For 
example, the stress acting on the interface may be important because it is often cou¬ 
pled to some interesting physical phenomenon, such as erosion, mechanotransduction, 
or movement of the entire porous blockage. However, in many models the behaviour 
of the flow near the interface is ignored, leading to incomplete information about the 
flow. As the lubrication equations used to approximate the flow away from the inter¬ 
face are reduced in order, approximated interf acial boundary conditions must be applied 
to couple the flows away from the interface (ICummings fe Waters! [20071 : 1 Waters et al. 


20061 : Cummings et al. 20091 : O’Dea et a,l\ 2009T 2013th A detailed understanding of the 


flow behaviour close to the interface will enable the correct coupling conditions to be 
applied and (dynamic) interfacial effects to be accurately included. It is, therefore, of 
fundamental interest to understand the flow near the interface of a porous blockage. 

Our main motivation arises from a tissue engineering problem. One approach to in 
vitro tissue engineering is to seed cells onto a porous biomaterial scaffold, which is then 
cultured within a bioreactor. The combination of cells and scaffold is referred to as a 
tissue construct and one particular aim of tissue engineering is to make porous tissue 
constructs as permeable as possible whilst maintaining their structural integrity. This 
is to enhance nutrient delivery to the cells residing in the interior of the construct via 
advection. Moreover, many cells are mechanosensitive, and their proliferation rate de¬ 
pends on the stress that they experience. In a high aspect ratio vessel bioreactor, the 
porous tissue construct is placed within a bioreactor shaped like a petri dish turned on 
its side, and the entire set up is saturated with a nutrient-rich fluid, and rotated around 
the bioreactor axis (figure [l]). The construct moves according to the force it experiences 
and, in the short term, the construct undergoes periodic motion. However, over a long 
time, the tissue construct drifts from its periodic orbit, an effect that can be attributed 
to weak inertia. To predict the construct trajectory and, ultimately, the nutrient trans¬ 
port in such a bioreactor, we must determine the forces acting on the construct. Thus, 
determining the flow through the porous construct is of particular importance. To gain 
insights into this problem, we consider the flow past a tissue construct held at a fixed po¬ 
sition, with a particular view to investigating the effect of weak inertia. Such an analysis 
can then be used to determine the construct trajectory when it is free to move. Another 
application of interest is the interfacial erosion of a porous biofilm, which is propo rtional 
to the square root of the fluid shear stress evaluated at the interface ( Duddu et al. 20091) . 
Thus, it is of interest to determine both the internal stress within the porous obstacle for 
tissue engineering applications and the interfacial stress acting on the obstacle boundary 
for biofilm erosion problems. 

The main aim of this paper is to investigate steady laminar flow through a narrow 
two-dimensional channel blocked by a tightly fitting finite-length porous obstacle, and 
then to generalise this analysis to investigate a time-dependent three-dimensional flow 
past a porous cylinder with an arbitrary smooth cross-section within a Hele-Shaw cell, as 
illustrated in figureEJ In particular, we are interested in the behaviour of the flow near the 
interface between the single-phase and porous regions, governed by the Navier-Stokes and 
Darcy equations, respectively. Transition from plug to Poiseuille flow in a two-d i men sional 


cha nnel is a classic problem, with the asymptotic structure exami ned by Van Dvkc (1970) 


Wilson dl971 ), and nume rical soluti ons given by , for ex ample, Brandt fe Gillis 


19661) 


and 
and 

would have to be used experimentally to induce uniform flow at the entry to a channel, 


Bodoia fc Osterlel ( 196ll ). Although IVan Dvkf ( 197fll ) mentions that a porous mesh 
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Figure 1. Schematic of a high aspect ratio vessel bioreactor, containing a porous scaffold (darker 
region). The bioreactor has a similar shape to a petri dish turned on its side, (a) Face view, (b) 
Side view showing the gaps between the scaffold and the bioreactor. Cells are seeded within the 
scaffold, and the entire bioreactor is filled with a nutrient-rich fluid and rotated around its axis, 
causing the scaffold to move. 


the nature of the coupled flow between single-phase and porous regions remains an open 
q uestio n - on e which we will answer in this paper. 

Thompson ( 19681) used the method of matched asymptotic expansions to analyse the 


three-dimensional flow near a solid circular cylinder within a Hele-Shaw cell, exploiting 
the small aspect ratio of the cell. Due to the no-flux condition on the solid cylinder, the 
normal flow near the interface is small. We are motivated by the tissue engineering appli¬ 
cation in which there is flow within a highly permeable porous obstruction of arbitrary 
smooth cross-section. The normal flow, in general, is then significant near the interface 
and the flows inside and outside the obstacle are coupled through their boundary condi¬ 
tions at the interface. Although we are interested in a highly permeable obstacle, we do 
not consider the Brinkman extension to the Darcy equations in this paper because the 
porous construct must maintain its structural integrity. That is, we consider a porous 
medium w hose underlying solid matrix is connected, so that Brinkman’s equations do 
not apply ( Auriault 20091) . 

We shall make extensive use of the method of matched asymptotic expansions by ex¬ 
ploiting the small aspect ratio of the channel/Hele-Shaw cell - one of our goals is to 
understand the flow near the interface (in an ‘inner’ region) given a general flow away 
from the interface (in an ‘outer’ region). Although the outer problems satisfy reduced 
equations, the inner problems are quasi-two-dimensional in each plane perpendicular to 
the interface. In particular, we determine the behaviour of inherently local properties, 
such as the stress acting on the interface. Additionally, we derive systematically the con¬ 
ditions that couple the outer equations for a general far-held forcing and for a general 
smooth cross-section of the porous obstacle. Due to their generality, the results obtained 
in this paper can significantly reduce the computational expense required to solve how 
problems with coupled single-phase and porous regions, whilst retaining important in¬ 
terfacial information. 

There is a large literature addressing the question of interfacial conditions on the 
boundar y of a p orous obstacle (and it remains an area of active research); see, for exam¬ 
ple, Nield fc Beian 12006 ) for a comprehensive review. We do not address this question 
here, and use as our interfacial boundary conditions: continuity of flux and continuity 
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Figure 2. The model geometry for flow within (a) a channel and (b) a Hele-Shaw cell. The 
light region is the exterior fluid and the dark region is the porous obstacle. 


of pressure, together with a no-tangential-slip condition. The first two are derived in 
Levy & Sanchez-Palencia (!1 975f ). the last is a special case of the general tangential slip 
condition derived in Carraro et all 1 20151) . 

The structure of this paper is as follows. In we consider the steady two-dimensional 
case of flow through a channel containing a tightly-fitting highly permeable porous obsta¬ 
cle. In 1 )2.11 we formulate the mathematical problem. In T2.31 we describe the asymptotic 
structure in the small-aspect-ratio limit and present the problems to be solved to couple 
the outer single-phase and porous regions. In l )2.4l and 1 )2.51 we investigate the asymptotic 
behaviour of these coupling conditions in the small- and large- Reynolds number limits, 
respectively. In particular, we find that the large-Reynolds-number limit induces a further 
boundary-layer structure, and investigate this in full. In lJ3]we extend our analysis to an 
unsteady three-dimensional flow in a Hele-Shaw cell past a tight-fitting highly permeable 
cylindrical porous obstacle. If the curvature of the cross-sectional boundary of the porous 
obstacle is not large (in a sense to be made precise later), the boundary layer structure 
in each plane perpendicular to the interface is equivalent to the two-dimensional case 
considered in and we are able to generalize the two-dimensional results to the three- 
dimensional case. In we apply the general results of to a time-dependent forced 
far-field flow in a Hele-Shaw cell containing a porous cylinder with circular cross-section. 
We determine the internal and interfacial stress within and on the porous construct, and 
show that fluid inertia can result in the generation of a net force on the cylinder, even 
when the far-field flow is periodic with zero mean. If the cylinder were allowed to move, 
this would cause a long-term drift, thus highlighting the singular perturbative nature of 
fluid inertia, and how it can be important in dynamical problems arising in lubrication 
flow. 


2. Steady two-dimensional flow 

2.1. Model formulation 

We consider the steady two-dimensional flow of an incompressible Newtonian fluid with 
density p and viscosity p in an infinitely long channel of height h. The channel contains 
a fully saturated porous obstacle with rectangular cross-section, height h and length 21, 
as illustrated in figure [21(a). The aspect ratio e = h/l is assumed to be small. We use 
Cartesian coordinates ( x , z) along and across the channel, as illustrated in figure^!, such 
that in the channel z € [0, h] and the plug lies in x £ [—1, l\, z £ [0, h\. The flow is driven 
by a prescribed upstream pressure gradient, such that in the far-field we have Poiseuille 
flow with cross-sectionally averaged velocity of magnitude t/ av . 

The flow exterior to the porous plug (the light regions in figurea)), which we refer to 
as the ‘exterior flow’, is governed by the incompressible steady Navier-Stokes equations 
with no body force. The exterior fluid velocity, pressure, and stress tensor are denoted 
by it = ue x + we z , p, and a, respectively, where e x and e z are the unit vectors in the 
x- and z-directions, respectively. The fluid inside the porous plug (the dark region in 
figure^), which we refer to as the ‘interior flow’, is governed by the incompressible 
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Darcy equations, again with no body force. The Darcy (volume-averaged) velocity and 
(interstitial) pressure within the porous plug are denoted by Q = Ue x + We z and P, 
respectively. The porous plug has constant permeability k. We do not use the porosity 
of the plug as a parameter as it is absorbed into the definition of the volume-averaged 
velocity, but note that it is assumed to be constant. 

We nondimensionalize by setting x = lx ', z = elz', u = U av u', U = U av U', w = eU av w', 
W = eU av W', p = {pU av /{eh))p' , P = (p,U av /eh)P', and a = (pU av /{eh))a' , where 
primes denote dimensionless quantities. Since all variables are dimensionless henceforth, 
we drop the primes without risk of confusion. 

The dimensionless governing equations in the exterior region |x| > 1, 0 < 2 < 1 are 
the scaled Navier-Stokes equations given by 

eRe (uu x + wu z ) = -p x + u zz + e 2 u xx , (2.1a) 

e 3 Re ( uw x + ww z ) = -p z + e 2 w zz + e 4 w xx , (2.16) 

0 = u x + w z , (2.1 c) 

where Re = pU av h/p is the modified Reynolds number based on the channel height. The 
governing equations in the porous plug |x| < 1, 0 < z < 1 are the scaled Darcy equations 
given by 

U = —KP X , e 2 W = —KP Z , 0 = U x + 114, (2.2a - c) 

where the dimensionless permeability K = k/h 2 . For the Darcy equations to be valid, 
we must have /{ < 1. We discuss the typical sizes of the three dimensionless parameters 
appearing in our model, e, Re , and K , in 1 12.21 


2.1.1. Boundary conditions 

On the channel walls we impose no-flux and no-slip conditions on the exterior flow and 
a no-flux condition on the Darcy flow, so that 


u = 0 on z = 0,1 for \x\ > 1; W = 0 on z = 0,1 for \x\ < 1. 


(2.3a, b) 


On the inflow (x = —1) and outflow (x — 1) inter faces, continuity of flux, conti¬ 
nuity of pressure, and a no-tangential-slip condition (ILevv fe Sanchez-Palencial 11975 ; 


Carraro et al. 20151 ) are given by 

u = U, p = P, w = 0 on x = ±1 for 0 < 2 < 1. 


(2.4a — c) 


The flow is driven by imposing a constant, upstream pressure gradient in the far-field, 
with 

dp 

-- > —12 as x —> — 00 , (2.5) 

dx 

ensuring that the cross-sectionally averaged horizontal component of velocity has mag¬ 
nitude 1 in the far-field. 


2.2. Parameter values 

Within a bioreactor, values of the three dimensionless parameters e, Re, and K can vary 
greatly in magnitude (table [T|. We will investigate the dimensionless problem defined 
by m (EH) in specific regions of this parameter space. 

The dimensionless problem is characterised by two lengthscales: the obstacle half- 
length, 1, and the channel height, e. We proceed by considering the case in which e <C 1, 
the channel height then being much smaller than the half-length of the obstacle. 

The value of K must be small enough for the underlying solid matrix to be connected, 
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Parameter Range 


Reference 


h 

l 

k 

U, 

P 

P 


2.5 x 10 -4 - 6 x 10 -3 m 
10“ 3 - 10 _1 m 
10“ 12 - 10“ 6 m 2 
3 x 10" 9 - 3 x 10 -2 ms -1 
10 3 kgm -3 

8.9 x 10 -4 kgm -1 s _1 


_Pazzano_e^_aZ i (2000); i Pathi _e£_aT (2005) 
Pierre^^e^^a^ (2008) 

Pazzcmo_et_aL' (2000); - Pathi_ef_a7 i (2005) 
Pierre_et_al^ (2008) 

Sim acck & Advani (1996) &ucosk^_et_a7 (2004) 
Nabovati_et_aZ i (2009) 
i ZhAo_ef_al i (2005); i Pierre_ef_al i (2008) 
Qiung_e£_aZ i (2007) 


e 2.5 x 10 -3 - 6 x 10° 

Re 10 -6 - 2 x 10 2 

K 2.7 x 10 -8 - 8.3 x 1(T 2 


Table 1. Typical dimensional parameter values and the corresponding values of the dimen¬ 
sionless parameters. The absolute upper bound of 1/12 ss 8.3 x 1CF 2 for K can be obtained by 
considering unobstructed Poiseuille flow through a channel. We use the density and viscosity 
values of water at room temperature for p and p, respectively. 


resu lting in Darcy’s equations governing the flow within the porous medium (Auriault 
l2009h . but large enough to enhance nutrient transfer via fluid advection to the cells 
within the plug. There is an absolute upper bound of AT < 1/12, and this can be seen by 
considering unobstructed Poiseuille flow through a channel. In reality, K <C 1, and we 
work in this limit. We see in the three-dimensional case, considered in that K = O(e) 
corresponds to an obstacle which is impermeable t o the flow a. t leading order. Hence, this 
limit reduces to the leading-order problem given in Thompson! ( 1968h . who considered the 
three-dimensional flow near a solid circular cylinder within a Hele-Slraw cell. We consider 
instead the limit in which e <C K <C 1, corresponding to the regime in which the plug 
is lon g, thin, and per meable at leading order in e. Examples of long, thin tissue s include 
skin ( Lei et alJlioTl h urothelial tissue ( Gabouev et al\ fibo.Qh . and the cornea ( Su et a,l J 


2003). Although we do assume that K <Cl throughout this paper, we effectively treat 


K = 0(1) as e —0 in our following asymptotic analysis for mathematical expediency. 

We must also define the relative scaling of Re. There is a distinguished limit when 
Re = 0(1), and the investigation of this limit is the subject of the next section. Further 
analytic progress can be made when Re is either small or large: we consider the asymptotic 
sub-limit Re <C 1 in 1)2.41 and the asymptotic sub-limit 1 <C 1/K <C Re <C 1/e in 1)2.5) 
We note that Darcy’s law (12.21) and the continuity of pressure condition (12.4b ) would 
require modifi cation if we were to generalize this work to significantly larger Reynolds 
numbers (Kayiany 120121 ). 


2.3. Asymptotic structure for Re = 0(1) 

We seek a solution using the method of matched asymptotic expansions with small pa¬ 
rameter e. The asymptotic structure is shown in figure[3ji, where each of seven asymptotic 
regions are labelled I-VII. We describe the regions as follows, noting that in each case 
the relevant transverse lengthscale is the channel height, e. There are three ‘outer’ re¬ 
gions characterised by an 0(1) axial lengthscale: two for the exterior fluid (regions I 
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Figure 3. The boundary-layer structure for Re = 0(1). The light region is the exterior fluid 
and the dark region is the porous obstacle, (a) The different asymptotic regions are indicated 
and labelled in Roman numerals, (b) A schematic representation of the flow. 



and VII) in which we have Poiseuille flow at leading order, and one for the interior fluid 
(region IV) in which we have plug flow at leading order. Thus, there is a transition from 
Poiseuille to plug to Poiseuille flow, between regions I, IV, and VII, respectively, as il¬ 
lustrated schematically in figure [3j>. The inflow transition from Poiseuille to plug flow 
occurs near the interface x = — 1, where there are two ‘inner’ regions, one on each side of 
the interface, each characterised by an axial lengthscale of 0(e). The exterior inflow inner 
region is region II, and the porous inflow inner region is region III. In a similar manner, 
the outflow transition from plug to Poiseuille flow occurs near the interface x = 1, where 
again there are two inner regions, one on each side of the interface, and each characterised 
by an axial lengthscale of 0(e). The porous outflow inner region is region V, and the 
exterior outflow inner region is region VI. 

2.3.1. Outer regions 

In the outer regions regions I, IV, and VII, we begin by posing asymptotic expansions 
of the form 

(u.Q,p,P) = (u 0 ,Q 0 ,p 0 ,P 0 ) + e(u 1 ,Q 1 ,p 1 ,P 1 ) + 0(e 2 ) as e-)• 0. (2.6) 

Substituting (12.61) into G3D, we find that the leading-order governing equations for the 
exterior flow are simply the lubrication equations 

0 = -pox + u 0zz , 0 = Pozi 0 = u O x + w O z for |x| > 1, 0 < z < 1. (2.7a, 6, c) 

Substituting (12.61) into m , the leading-order governing equations for the interior flow 
are 

U 0 = -KP 0x , 0 = -KP Oz , 0 = U 0x + W 0z for |x| < 1, 0 < z < 1. (2.8a, b, c) 

At leading-order, the channel wall boundary conditions (12.31) become 

Uq = 0 on z = 0,1 for |x| > 1, (2.9a) 

Wo = 0 on 2 = 0,1 for \x\ < 1, (2.9 b) 

while the far-held condition (12.51) becomes 

—> —12 as x —> — oo. (2-10) 

dx 

We are unable to impose at leading-order the interfacial boundary conditions (12.41) in the 
outer regions; these conditions are satisfied within the inner regions described in 1 )2.3.21 
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For the outer problems, we mu st instead impose matching conditions (using the matching 
law given in i Van Dvke (1 9751 )). These are determined in 92.3.21 and we state them here 
to close the leading-order outer problem. As the outer pressure is independent of z and 
the fluid velocity is proportional to the pressure gradient, we find that we may impose 
continuity of total flux and pressure at each interface (corresponding to Neumann and 
Dirichlet conditions on the pressure, respectively). That is, 


Uq dz = / Uq dz, Po = Pq on x = ±1. 


(2.11a, b) 


We note that whilst the conditions (12.111) close the outer problem, they give no informa¬ 
tion about, for example, the leading-order stress acting on the interfaces. 

The leading-order solutions are readily obtained as follows: in region I, 


«0 = 

6z(l - 

- z)e x , 

Po = 

— 12a;; 

(2.12a, b) 

in region IV, 






u 0 - 

— &X ? 

Po = - 

1 + X 

K 

+ 12; 

(2.13a, b) 

and in region VII, 






CO 

II 

O 

3 

z)e x , 

Po = 

— 12a; -f 

-Kiy 

(2.14a, b) 


The pressure is determined up to an arbitrary constant which (at this order) we fix by 
supposing p 0 + 12a; —» 0 as a; —>■ — oo, without loss of generality. We remark that the 
presence of the porous plug introduces a jump in pressure of 24 — 2 /K ~ —2 /K (as K 
is small). Equivalently, an 0(1) pressure drop creates an 0(K) flow through the porous 
p lug. 

We note that the horizontal components of velocity (12.12h ) and (I2.13h ) are incompat¬ 
ible with the continuity of flux condition (12.4k ') on the interfacial inflow boundary for 
the full problem. Similarly for outflow, we find that (12.13b l and (12.14b ) are incompatible 
with (12.4b .). In the next section we introduce inner regions to resolve this issue, but for 
the remainder of this section we consider higher-order terms in the outer regions as the 
higher-order pressure is coupled to the leading-order flow in the inner regions. 

Using the leading-order solutions (12.121) and (12.141 ) . the 0(e) terms in the exterior-flow 
equations (EH) become 

o = —pix + Uizz, 0 = Pizi o = Mia, + Wiz for |a;| > 1, 0 < z < 1. (2.15a-c) 

Similarly, using the leading-order solutions (12.131) . the 0(e) terms in the interior-flow 
equations (12.21) become 

U 1 = -KP lx , 0 = -KP lz , 0 = U lx + W lz for \x\ < 1, 0 < * < 1. (2.16a - c) 

The channel wall boundary conditions (12.31) at this order are as follows 


Mi = 0 on 2 = 0,1 for |x| > 1, 

(2.17a) 

Wi =0 on 2 = 0,1 for \x\ < 1, 

(2.176) 

while the far-held condition (12.51) implies 


dpi n 

—- > 0 as x —> —oo. 

da; 

(2.18) 


The matching conditions across each interface are determined in 92.3.21 and we give 
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them here to close the problem. For inflow, we have 


f u±dz= f Uidz, pi — Pi = II (Re,K) on x = — 1, 
Jo Jo 

and for outflow we have 

f Uid 2 = f Uidz, pi — Pi = II + (i?e, K) ona) = l, 
Jo Jo 


(2.19a, b) 


(2.20a, b) 


where the pressure jumps II - and II + are functions of the parameters Re and K. In 
1 )2.3.21 we specify the problems that need to be solved to determine these functions, 
while in 1)2.41 and 1)2.51 we determine their asymptotic behaviour for small and large Re. 
Given II^, the solution to (12.151) (12.201) is readily determined: in region I, 


in region IV, 


and in region VII, 


Ui = 0, 

Pi = n ; 

(2.21a, b) 

Qi = o, 

Pi =0; 

(2.22a, b) 

Ml = 0, 

pi = n+. 

(2.23a, b) 


As with the leading-order solutions (12.121) (12.141) . the pressure field is unique up to an 
arbitrary constant. For convenience, we choose the constant such that Pi(0,0.5) = 0, 
without loss of generality. We proceed by considering the inner regions, highlighting the 
problems that would need to be solved numerically to determine the functions II - and 

n+. 


2.3.2. Inner regions 

The appropriate scalings in the inflow and outflow inner regions are x = 4=1 + eX^, 
respectively, with (w, W ) = ( w , W)/e. To make it clear that we are working in the inner 
regions we introduce overbars on the other dimensionless variables, writing ( u , U ) = 
(u, [/), (p,P) = (p, P) and (m,Q) = ( u,U)e x + ( w,W)e z . Using vector operators in 
terms of and z as appropriate, the exterior fluid governing equations ( 12 . 11 ) become 


Re (u • V) u = — e 1 Vp + V 2 tt, 0 = V • u in regions II and VI, (2.24a, b) 

which are valid when X~ < 0 for inflow (region II) and X + > 0 for outflow (region VI), 
with 0 < 2 < 1 in both cases. The interior fluid governing equations (12.21) become 

Q = —e _1 AVP, 0 = V • Q in regions III and V, (2.25a, b) 

which are valid for X~ > 0 for inflow (region III) and X + < 0 for outflow (region V), 
with 0 < ^ < 1 in both cases. 

Similarly to (12.fill . we form inner expansions in powers of e as follows: 

(u,Q,p,P) = (uo,Qo,Po,Po) + e(ui,Qi,p 1: Pi) + 0(e 2 ) as e ->• 0. (2.26) 

At leading order, it follows from (12. 24b , ~) that 

0 = —Vp 0 in regions II and VI, (2.27) 

while (I2.25k l gives, similarly, 

0 = — KXPq in regions III and V. (2.28) 


The appropriate interfacial conditions on X ^ = 0 for the leading-order pressure are 
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Region II Region III z = 1 


uq = 0 , wq = 0 


M 0 ~ 6z(l - z)e x 
pi - 12X~ + IR 



Wq = 0 


Q o 


P i ~- 


X- 

~K 


Re (uo • V) u 0 = —Xp 1 + V 2 «o 

1 

1 

1 

1 

Q 0 = -KVP 1 

0 = V • U 0 

■ 

1 

1 

0 = V • Qq 


uq = 0, w o=0 


Wq = 0 


X- = 0 2 = 0 

Figure 4. The inner inflow problem for Re = 0(1). The flow is from left to right. 


given by the leading-order terms in (12.4b b which become 

Po = Pa on X T = 0. 


(2.29) 


Thus po and Pq are constant and equal in regions II and III and in regions V and VI, 
leading to the matching condition (I2.11b h Using the leading-order outer solutions (12.121) 
(12.141) for the pressure, we find 


Po = 12 in region II, Pq = 12 in region III, 


for inflow, and 


Po = 12 —— in region V, 
A 


p 0 = 12 — — in region VI. 


(2.30a, b) 


(2.31a, b) 


for outflow. We no te that these leadi ng -order pressures are constant on the interface, as 
predicted by iLevv fc Sanchez-Palenciaf ( 1975l l. 

As the leading-order fluid pressures are constant within each inner region, the leading- 
order inner flow is driven by the first correction to pressure, and these equations are given 
by the 0(1) terms in (12.241) and (12.251) . We illustrate the leading-order flow problems in 
figure [4] for inflow, and in figure [5] for outflow. In both cases, the inner flow transition 
problem is governed by the two-dimensional steady Navier-Stokes equations coupled to 
the Darcy equations. To fully determine the inflow and outflow pressure jumps II - and 
n+ , we would have to solve the problems indicated in figures |4] and [5] numerically for two 
parameters: Re and K. 

However, to couple the outer regions, we require another condition on the pressure. 
We note that we can analytically determine suitable coupling conditions for the hori¬ 
zontal components of velocity averaged over the channel height by appealing to global 
conservation of mass, rather than solving the full problem. By integrating the continuity 
equations across the channel height, using the conditions for the flow in the normal direc¬ 
tion at each boundary, and matching with the outer solutions, we deduce the following 
coupling conditions for the outer horizontal components of velocity 


/ u 0 dz = Uq&z, 
la Jo 


/ Ui dz — I U\ dz 

/ o Jo 


ons = ±l. 


(2.32a, b) 


We proceed by determining the inner flow in the asymptotic sub-limits Re <C 1 and 
1 < 1/K <C Re <C 1/e, starting with the former. This will allow us to calculate quantities 
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Region V Region VI z = 1 


Q o 

Pi - 


W 0 = 0 


a+ 

If 


uo = Uo 
Pi = -Pi 
w 0 = 0, 


W 0 = 0 


uq = 0, w o = 0 


Pi 


uo ~ 62(1 - z)e x 

-12A+ + n+ 


Q 0 = -KXP 1 

1 

1 

1 

1 

Re (uq • V) u Q = — Vp 1 + V 2 «o 

0 = V • Q 0 

■ 

1 

1 

0 

II 

<1 

SI 

c 


uo = 0, wq = 0 


X+ = 0 - = 0 

Figure 5. The inner outflow problem for Re = 0(1). The flow is from left to right. 


that are unobtainable from sole consideration of the outer problem, such as the stress 
acting on the interface and the limiting behaviour of the functions II~, II + . 

2.4. Small Reynolds number: Re -C 1 

In the sub-limit in which Re <C 1, the leading-order problems outlined in figures [4] and [5] 
become Stokes flow coupled with Darcy flow. Therefore, the flow equations are now 
reversible, and it suffices to solely consider inflow. We therefore use A'” = X for ease of 
notation. Further, we note that the inflow and outflow pressure jumps are now linked via 
the expression II + = —II - . 

We solve the resulting leading-order system numerically (for the streamfunction within 
the exterior flow and for the pressure within the porous flow) using a second-order accu¬ 
rate central finite-difference scheme. Due to the elliptic nature of the governing equations, 
we use an iterative method to bypass some of the cost of solving the fully coupled elliptic 
problems in one step. That is, we iterate between solving Laplace’s equation for the pres¬ 
sure in region III, and the biharmonic equation for the streamfunction in region II, and 
use the solution found in the previous iteration to determine the boundary conditions 
for the current iteration. We truncate the infinite domain to [—1,0] for X in region II 
and [0,1] for A' in region III (we note that an extended domain adds little accuracy to 
the scheme), using 300 grid points for A in each region and 300 grid points for 2 . The 
simulation is halted once the difference between successive iterations at each grid point 
is less than 10~ 4 . 

The normal and shear components of stress acting on the interface are 

Oxx ~ -12 + e(-pi +2u 0 x), (2.33 a) 

o X z ~ e (uoz + wox), (2.336) 

respectively, as e, Re —> 0. We show the O(e) stress results for K = 0.01 in figure [(3 
and calculate II” = 1.92 to three significant figures. The leading-order normal stress 
term in (12.33al) takes the specific value of —12 because we have imposed P(0, z) = 0 
to uniquely define the pressure. In fact, it is the O(e) terms in (12.3311 which are more 
interesting, as these are the dominant terms which vary in 2 and, unlike the leading- 
order pressure, cannot simply be determined from the outer solution. The boundary 
layer analysis is essential to obtain these results, which cannot be determined solely from 
the outer solutions. The reason the outer solutions cannot make good predictions of the 
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Figure 6. The 0(e) normal and shear stresses acting on the interface X = 0, for K = 0.01. 

(a) pi - 2uqx (b) - (m 0 z + viox)- 

^-variation of interfacial stress is because of the tendency of the exterior flow to move 
towards the channel edges in transition between Poiseuille and plug flow in the boundary 
region. The interfacial normal stress is greater/smaller than predicted by the outer flow 
towards the centre/edges of the channel, whilst the magnitude of interfacial shear stress 
is greater than predicted by the outer solution. 

2.5. Large Reynolds number: 1 -C l/if<C Re <C 1/e 

We now consider the sub-limit in which 1 <C 1/K <C Re <C 1/e. We do not consider 
further the sub-limits in which ReK = 0(1) or 1 -C Re -C 1/K -C 1/e. This is because 
we are interested in the tissue engineering problem with porous plugs whose values of 
K are not too small. Moreover, the above sub-limits yield coupled nonlinear problems 
which are not amenable to analytic study. We emphasize that the limit we consider is 
when the external fluid has a large Reynolds number, but the pore Reynolds number is 
small, so that the Darcy equations and continuity of pressure conditions are valid. 

In this section we perform an asymptotic analysis using the small parameter 6 = 
1/Re -C 1 for ease of notation. Therefore, when we refer to, for example, ‘leading-order’ 
in this section, it is meant in reference to 6. In the limit S —> 0, there is a significant 
change in boundary layer structure within the inner regions II and VI compared to the 
0(1) Reynolds number case. The limit as <5 —> 0 is singular, and we expect viscous effects 
to become important in regions near each boundary. Furthermore, as fluid momentum 
dominates within regions II and VI, the flow direction is important. Indeed, we obtain a 
different boundary layer structure for inflow and outflow, and therefore split the analysis 
into these two cases, starting with inflow. 

In the limit of small 5, the first-correction terms Mi, Qi,p\, and P± from the asymptotic 
series in the outer regions m are now scaled with S 1 to be consistent with the new 
inner scalings, as follows: 

(ui,Qi,Pi,Pi) = <5 -1 («io,<3io>Tio,-Pio) + (uu,Q n ,pu,Pii) + O(S). (2.34) 

The scalings for the pressure jumps II - and II + will turn out to be 

IT =11/+ 0(5), 11+= cT 1 ^ + 0(1) as5->-0, (2.35a, 6) 

our objective being to determine 11/ and Ilg". 

Though we have not yet fully described the inner asymptotic structure, we give the 
form of the asymptotic expansions within each region in Table [2] for ease of reference. 
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Region 

Asymptotic expansions 

Relative scalings 

I 

(u,p) = (u 0 ,po) + e(«i,pi) + Of 2 ) 

(ui,pi) = 8~ 4 (uio,pio) + (ttii,pii) + 0(5) 


II 

(u,p) = (u 0 ,po) + e(«i,pi) + Of 2 ) 
u 0 = 6z(l — z)e x + 5u 0 i + 0(S 4/3 ) 

Pi = Pu + 0(<5 1/3 ) 

u = u, w = e _1 ui 
x = eX, p = p 

Ha 

(u,p) = (uo,Po) + e{ui,pi) + Of 2 ) 
uo = 62(1 — z) + 5uoi(0, z) + 0(5 4/3 ) 
wo = woi + 0(5 1/3 ) 

Pi = Pn(0, z) + 0(5 1/3 ) 

u = u, w = e~ 1 8w 
x = e8X, p = p 

III 

(Q, p) = (Q 0 , Po) + e (Qi, Pi) + Of 2 ) 

Qo = Qao + 5Qoi + 0(8 4 ' 3 ) 

Pi = Pu + 0(5 1/3 ) 

U = U,W = e^W 
x = eX,P = P 

IV 

(Q, P) = (Q 0 , Po) + e (Qi, Pi) + Of 2 ) 
(Qi,Pi) = 5- 1 (Q 10 ,Pio) + (Qu,Pii)+0(8) 


V 

(Q, P) = (Qo, Po) + e ( Qi,Pi ) + Of 2 ) 

Q 0 ~ e x + 5 1 ' 2 Q 0a 

Pi ~ -X~/K + 5 1/2 Pi b 

h 

Itt, 

II II 

11 11 

H 

VI 

(u,p) = (uo,po) +e(ui,pi) + Of 2 ) 
u 0 ~ e x + 8 1/2 u 0 a 

pi ~ 8~ 1/2 p la 

u = u, w = t 4 w 
x = eX, p = p 

Via 

(u,p) = (uo,Po) + e(ui,pi) + Of 2 ) 
uo = uoo + 0 ( 8 1/2 ) 
wo = mi + 0(8 1/2 ) 
pi = <5 -1 pio + 0(8~ 1/2 ) 

u = u, w = e~ 1 8w 
x = e5~ 4 X, p = p 

VII 

(u,p) = (u 0 ,po) + e (ui,pi) + Of 2 ) 

(ui,pi) = 5- 1 (uio,pi 0 ) + 0(8~ 1/2 ) 



Table 2. Asymptotic expansions within boundary layers in the large-Reynolds-number limit. 


2.5.1. Inflow 

The inflow inner problem (within regions II and III) is shown in figure [U The limit 
as 5 —> 0 yields leading-order inviscid flow in region II away from the boundaries, with 
five additional boundary layers being required to satisfy all of the boundary conditions. 
Two of these viscous boundary layers are on the channel walls, one is on the porous 
interface, whilst the remaining boundary layers occur in the corners between the channel 
wall and porous interface. The boundary layer structure for the inflow problem is shown 
in figure [7] 
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0(e6 1 ' 3 ) 


0(e6 1 > 3 ) 


0(e) 


0(eS) 

—y 



Figure 7. Schematic diagram of the inner regions Ila-c within the inflow boundary region II. 
The light region is the exterior fluid and the dark region is the porous obstacle. The size of the 
boundary layers has been exaggerated for illustrative purposes 


We start by considering the coupled regions II and III, and show that the transition 
to Poiseuille flow occurs solely within region III at leading order in S. 

2.5.2. Regions II and III 

Within regions II and III, we make the following asymptotic expansions 

uq = 6z(l — z)e x + 5u 0 i + (5 4/3 m 02 + 0(6 5 ^ 3 ), (2.36a) 

Qo = Qoo + ^Qoi + S 4/3 Q 02 + 0(<5 5 / 3 ), (2.366) 

C Pi, Pi) = (Pn,Pn) + S 1/3 (p 12 ,Pi 2 ) + 0(S 2 / 3 ) (2.36c) 

as 6 —> 0 , where the non-integer powers of S come from corrections to the flow that arise 
due to the presence of a boundary layer near the channel walls (regions lib in figure [7]). 
These terms are used in Appendix [A] but are not required further in this section. 

Substituting the asymptotic expansions (12.361) into the governing equations presented 
in figure^ the 0(1/6) terms in region II are trivially satisfied, but so far Q 0 o is unknown. 
Proceeding with the expansions, we find at 0(1) that Uqi = (moi^oi) is governed by 

uooUqix + «oo z W()i = -Pnx ~ 12 > u 0 oW 0 ix = ~Pii z , V • u 0 i = 0 in region II, 

_ _ (2.37a-c) 

where u oo = 6z(l — z), while Q 00 = (Uoo, Woo) is governed by 

Q 00 = —K'VPii, 0 = V • Q 00 in region III. (2.38a, 6 ) 

The boundary conditions we impose at this order in region II are again obtained via 
matching. On the channel walls (via matching with region lib in Appendix [X]) , we have 
no flux, so that 

wq\ = 0 on z = 0,1 for X~ < 0, (2.39) 

while in the far-held we match with region I, to obtain 

Moi —y 0, Jin ~ —12A' - + Ilf as X~ —> —oo for 0 < z < 1. 


(2.40a, 6) 
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Similarly, the boundary conditions in region III are 

W 00 = 0 on z = 0,1 for X~ > 0 , 
while matching with region IV gives 

_ 

Pn ~ —as X~ —> oo for 0 < 2 < 1 . 


(2.41) 


(2.42) 


We note that the far-field matching condition Q 0 o ~ e-x as X~ —» oo follows consistently 
from (12.4211 using (12.38b ,). 

The interfacial conditions at X~ =0 for (12.371) (|2.38D are formally derived from con¬ 
sideration of the inner-inner boundary layer (region Ha) in 1 12.5.31 within which we are 
able to satisfy the no-tangential-slip condition, but additionally find that the leading- 
order flow and the first correction to the pressure pn are unchanged. From 1 12.5.31 the 
correct conditions to couple regions II and III are continuity of flux and pressure, with 

6z(l — z) = Uqo, Pn = P\i on X~ = 0 for 0 < z < 1 . (2.43a, b) 


The system to solve in regions II and III (up to the 0(6) terms) is given by (12.371) 
(12.431) . We see that regions II and III have decoupled: the boundary condition (|2.43b ) 
allows us to solve for Q 0 o and Pn using (|2.38D with (12.411) (12.4211 . and then we can use 
the boundary condition (12.43b ") to determine u 01 and pii using (12.371) with (12.391) (12.401) . 
The governing equations (12,38f) . with boundary conditions (12.411) (12.421) and (12.43b '). are 
solved by 

°o „ 

—KPn = X~ + ^ Ak exp(—2fc7rV~) cos( 2 £; 7 r,z), Ak = - 3 . (2.44) 

fc=1 k ) 


Since Uqo is independent of X~ within region II, the leading-order transition from 
Poiseuille to plug flow occurs entirely within the porous medium in region III, and is 
governed by (12.441) . 

We now complete the solution in region II by determining u 01 and pn. As well as 
allowing us to determine IIJ”, knowledge of pn gives the normal stress acting on the 
interface up to O(e), the first order at which the stress varies in the ^-direction. A 
rearrangement of the exterior flow equations (12.371) allows us to decouple the region II 
system (where — 00 < X~ < 0 and 0 < z < 1) into 

z( 1 — z)X 2 woi + 2woi = 0, (2.45a) 

V 2 pn = -2u 0 ozWoix- ■ (2.456) 

We solve (I2.45al) for Wq 1 first, using boundary conditions (12.391) on z = 0, 1 and 
(12.40b ,') as X~ — 00 . We derive a consistent boundary condition for uioi on X~ = 0 
using the boundary condition (12.43b ) for the pressure on X~ = 0, in combination with 
the governing equation (12.37b ). giving 


w oix - 


1 

Kuoo 


OO 

2knAk sm(2kTrz) 
k =1 


on A =0 for 0 < z < 1 . 


(2.46) 


The boundary condition (12.461) ensures that Wqi does not vanish on X~ = 0, induc¬ 
ing a slip velocity on this interface within region II. This issue is taken care of in an 
inner-inner boundary layer closer to the interface (region Ha in figure [7J in 1 12.5.31 The 
linearity of the governing equation (12.45 a|) . with boundary conditions (12.391) . (12.40^ . and 
(12.461) . also allows us to deduce that woi scales with K -1 . Further, from the continuity 
equation (12.37b ) it follows that uqi also scales with K . 
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(b)i 



0.5 


-0.5 

X~ 


Figure 8. The streamlines for inflow near the interface, (a) The streamlines for the scaled 
first-correction velocity Kuoi, which is independent of K. The fluid moves from the centre of 
the channel to the channel walls, (b) The streamlines for the full velocity up to first-correction 
uoo + 5uoi (with a large value of 5/K = 1 taken to emphasise the effect of the correction). The 
fluid moves from left to right. 


We solve for w oi using a standard central finite-difference scheme. By numerically 
integrating woi with respect to A' - , we are able to plot the streamlines in figure [8] 
In figure [8b, we show the streamlines for the first-correction to the velocity, and in 
figure [8b, we show the streamlines for the full velocity up to this first-correction. We 
see that the fluid moves from the centre of the channel towards the channel walls, as 
expected. Therefore, adjustment from Poiseuille to plug flow occurs in region II, but 
at higher order. This movement induces a horizontal slip velocity on the channel walls, 
which is resolved in Appendix [A] by introducing a boundary layer near the channel walls 
in which z = 0(A 1 / 3 ) on the bottom wall, with a similar layer near the top wall. These 
regions are labelled as lib in figure [7] In Appendix [A] we show that the flow near z = 0 
has the form 

_ , Ai(afcv) du , 

C k exp(p k X )g k { 0) Jo oo - as z->• 0+, (2.47) 

J 0 Ai(o fe s) ds 

where Ck, Hk , and gk are described in (lA~T1) (1A21) . Ai is the Airy function of the first 
kind, and a*, is defined in (I A 10b ). From (12.471) . we see that the leading-order solution is 
preserved near the channel wall. In fact, the slip is resolved within higher-order terms in 
the boundary layer, and does not affect the leading-order solution in region II. 

The problem for pn is given by (12.45 6[1 with boundary conditions (12.40b ) and (|2.43b ). 
We obtain consistent boundary conditions for pn on the channel walls by substituting 
the leading-order flow solution into the governing equation (12.37b ) to obtain 

p llz = 0 on 0 = 0,1 for X~ < 0. (2.48) 

We solve the full system for pn, given by (12.4561) with boundary conditions (12.40b ). 
(12.43b ). and (|2.48[) numerically using a central finite-difference scheme. Using a similar 
scaling argument as for m)oi, we find that II)” scales with 1/K. We are able to calcu¬ 
late II j” by plotting (pn + 12A - ) K, and using the matching condition (12.40b ). As IIj” 
is independent of z, and to more easily show the pressure drop on a graph, we plot 


Mo ~ 6z(l — z) + 6 

fc=l 
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determine the outer pressure drop for inflow. 


Pn dz + 12X j K for varying X in figure [9] We find thereby that n x = a/K, 


where a = —0.117 to three significant figures. 


2.5.3. Inner-inner region Ha 

We now investigate the inner-inner boundary layer closer to the interface, region Ha 
in figure [3 This boundary layer is required to satisfy the no-tangential-slip condition on 
the interface, and captures a balance between the relevant inertial terms and the viscous 
terms induced by the presence of the porous interface. 

The relevant scalings are X~ = 6X and Wq = Swo, and we introduce the new variables 
uq = uq and pi = p\ to signify that we are working in region Ha. Under these scalings, 
in region Ha (where —oo < X < 0 and 0 < z < 1) the exterior flow equations in figure 0] 
become 


S SwqUqz — PlX d - ^ Y (2.49 a) 

uow 0 x + S 2 w 0 w 0z = -5pi z + w 0 xx + fi 2 ™0zz, (2.496) 

0 = u 0 x + S 2 w 0z , (2.49c) 

and the interfacial boundary conditions become 

uq = Uo, pi = Pi, wo = 0 on X = 0 for 0 < z < 1. (2.50a — c) 

The purpose of this boundary layer is to determine the tangential velocity wo, and to 
ensure that it is possible to satisfy the no-tangential-slip boundary condition on X~ = 0. 
Therefore, we form asymptotic expansions as follows: 

Mo = 6z(l — z) + O(S), (2.51a) 

wo = w 0 i + 0(S 4/3 ), (2.516) 

pi =Pn(0,z) + 0(S 1/3 ) (2.51c) 


as 6 —> 0, where p n is known from 1)2.5.21 Substituting the asymptotic expansions (12.511) 
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into the equations (12.4911 and equating powers of <5, we find that (12.49a.li and (12.49d) are 
automatically satisfied at leading order, whilst the leading-order terms in (12.49&I) are 

6z(l — z)w Q1 x = ® 01 jx' (2.52) 

Using (12.441) . the leading-order interfacial boundary conditions (12.501) become 

6 z(l — 2 ) = —KP 1X y, p lx = Pn, u?oi =0 on X = 0 for 0 < 2 < 1. (2.53a — c) 

The relevant far-field condition arises via matching with region II, giving 

iuoi —> Tuoi(0,z) as X —> —00 for 0 < 2 < 1, (2-54) 

where wqi is determined in 112.5.21 

The interfacial conditions (12.53h .bl are used to couple regions II and III in 112.5.21 The 
remaining terms in the system (12.521) (I2.54|) yield w 01 , and are solved by 


w 0 i = ui O i( 0 , 2 ) (l - exp( 62 (l - z)X)J 


(2.55) 


We see that (12.551) can only be a solution for inflow, as we require 62(1 — z)X < 0 as 
we move back into region II. For outflow, this inner-inner region is not a valid boundary 
layer, and we must tackle the problem in a different manner. 

We note that this boundary layer persists for alternate tangential slip conditions. We 
can showcase the robust nature of the boundary layer vi a two ex am ples . The first example 
uses the general tangential-slip condition derived in Carraro et all ( 20151) . That is, instead 
of (12.53b ). we consider a condition of the form 


w = Xu on X = 0, for 0 < 2 < 1, 


(2.56) 


which reduces to (12.53b ) when A = 0. Here , A is a slip coefficient, obtained by solving the 
cell problem given in lcarraro et al. ( 2015?) . Using this condition, the transverse velocity 


is 


w ~ 6 A 2(1 — 2 ) exp( 62 (l — z)X) + Ju>oi(0, 2 ) ^1 — exp( 02 (l — 2 )X)j + 0(5 4 ^ 3 ), (2.57) 

in region Ha. In the second example, we use a Beavers and Joseph slip condition dBeavers fc Joseph 
19671) of the form 


A 

y/K 


(W — w) on X =0, for 0 < 2 < 1, 


(2.58) 


where A is the experimentally determined slip coefficient. In this case, the transverse 
velocity would be 


m ~ J (®oi(0, z) + 6 ^°_ 2) ) + W 4/3 ). <2.5S» 

Hence the boundary layer structure persists if the no-slip condition is modified in either 
of these two ways. Physically, we can conclude that the thickness of each boundary layer 
involved in inflow is unchanged for any of the three boundary conditions applied at the 
interface. Moreover, the flow induced by each of the three boundary conditions is only 
different in region Ha, the boundary layer closest to the interface. 

Returning to the original no-slip boundary condition (12.50b ). it is straightforward to 
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Figure 10. The interfacial stress predicted by the inner flows (from (12.6011 1 at 0(e). (a) The 
normal stress (we plot Kp n(0, z )) and (b) the shear stress (we plot 6z(l — z)Kwoi (0, z)). As both 
the normal and shear stress scale with 1/K at this order, both plotted terms are independent 
of K. 

calculate that the components of stress on the interface at x = — 1 are given by 

Oxx = - (12 + ep n (0, z)) + 0(S 1/3 e), (2.60a) 

u xz = -e (6z(l - z)wqi (0, z)) + 0(S 1/3 e) (2.606) 

as e,<5 —»• 0, with e <C S <C 1. Whilst the normal stress (12.60al) is nominally 0(1), the 
leading-order term comes from the decision to impose P(0, 0.5) = 0 to uniquely define 
the pressure. Thus, the more interesting result occurs at O(e), which is the lowest order 
at which the normal stress varies in the z-direction. For the shear stress (12. 60611 . the 
leading-order terms are O(e), and this varies in z. Those terms which vary in z cannot 
be determined solely from the outer solution: they are obtained from our boundary layer 
analysis and are each plotted in figure [10] We show in the previous section that pn and 
woi both scale with 1/if, thus the stress also scales with 1/if at O(e). 

Additionally, since we have imposed P(0,z) = 0, we have no 0(e/S) constant term 
from the pressure in the normal stress (I2.60al) . In three-dimensions, the tangential (that 
is, in the (e x x e z )-direction) variation in pressure ensures that there will be an 0(e/5) 
term in the corresponding normal stress component no matter how we define the pressure 
(as we show in £ 13.6.11) . 

We have now solved for the flow variables within the inflow boundary layers for the re¬ 
gions which affect the stress and the pressure drop. As shown in figure [3 there are further 
boundary layers which are important for inflow, and we discuss these in Appendix El 

2.5.4. Summary 

The main part of the leading-order inflow may be summarized as follows. The leading- 
order axial flow and pressure are unchanged until they reach the interface, after which 
they transition to plug flow inside the obstacle in region III, whose length is comparable 
to the channel height. This has an effect on the correction to normal flow and pressure 
in region II. Meanwhile, the tangential velocity satisfies the no-tangential-slip condition 
on the obstacle surface in region Ha - a boundary layer of even smaller width, which 
is comparable to the length of the channel height multiplied by the reciprocal of the 
Reynolds number. 

Combining the original asymptotic expansion (12.261) with the asymptotic expansion in 
region II, namely (I2.36d) . we deduce that the pressure jump between outer inflow regions 


















20 M. P. Dalwadi, S. J. Chapman, S. L. Waters, and J. M. Oliver 

I and IV is 0(e). Specifically, we have determined that 

IT ~ as e, <5 —> 0 with e<i5cl, (2-61) 

where a ~ —0.117. 

The argument presented in this section will not hold if the fluid is leaving the obstacle 
rather than entering it. This is because the boundary layer region Ha is only valid for 
inflow (as can be seen from (12.55b f ). We reconsider our approach for outflow in the next 
section. 

2.5.5. Outflow 

The outflow inner problem (within regions V and VI) is shown in figure [5] We again 
seek a solution as <5 —► 0 with e <C S -C 1 using the method of matched asymptotic 
expansions. 

The asymptotic structure for outflow is shown schematically in figure 1111 The flow 
in region V is unchanged at leading order and reaches the interface as uniform plug 
flow with magnitude 1. There is an inviscid core in the centre of region VI, flanked by 
Prandtl boundary layers near z = 0,1 whose thicknesses are of 0((<5V + ) ' ). These 
boundary layers (region VIb) grow downstream, and the transition to Poiseuille flow 
occurs when their thicknesses become of 0(1). This occurs on a horizontal lengthscale 
of 0(1/8) in a transition region we denote region Via. There is a third boundary layer, 
denoted region Vic, which occurs in the corner between the channel wall and the porous 
interface. This region deals with a singularity at the corner that arises in region VIb, as 
in the classic Prandtl problem. 

In this section, we investigate the regions VI and Via to understand the transition from 
plug to Poiseuille flow and to calculate the pressure jump between the outer regions IV 
and VII. Regions VIb and Vic are discussed in Appendix IB.ll since they are not important 
to the transition. Recall that, from (12.611) . we calculated an O(e) pressure jump between 
outer regions I and IV for inflow. We show that the pressure jump between outer regions 
IV and VII is of 0(e/8) for outflow. 

We start by investigating how the fluid moves through region VI. The leading-order 
behaviour is plug flow, and we determine the effect of the interfacial no-tangential-slip 
condition on the correction to the plug flow. In particular, we determine the far-held 
behaviour of the flow, which allows us to impose the correct matching conditions for the 
problem of transition to Poiseuille how in region Via. 

2.5.6. Inner regions V and VI 

To determine the outhow pressure jump (I2.35b h we must consider hrst-order how 
terms. Therefore, we pose the asymptotic expansions 

uo ~ e x + S 1 / 2 u Qa , Qq - e x + S 1/2 Q 0a , p 1 ~ S~ 1 / 2 p la , Pi -^ + 8 1 / 2 P lb , 

(2.62a- d) 

as 8 —> 0, where the non-integer powers of 8 arise due to the correction terms from the 
boundary layers in region VIb near the channel walls. We determine the leading-order 
behaviour within the coupled regions V and VI, but do not fully solve for the higher-order 
terms due to their complexity. Instead, we analyse the higher-order problems to obtain 
the information required for the matching condition with the transition region Via. 

Substituting the asymptotic expansions (12.621) into the system outlined in figure O the 
leading-order equations are trivially satisfied. The 0(8 1 ^ 2 ) exterior flow equations are 

UoaX+ = ~PiaX+ j WoaX+ = ~Pi a z, 0 = u 0aX + + WQaz in region VI, (2.63a - c) 
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Figure 11 . Schematic diagram of the inner regions for outflow - Via, VIb, and Vic. The size of 
the boundary region VIb has been exaggerated for illustrative purposes, is of height 0(e<5 1//2 ), 
and grows proportional to the square root of the distance from the porous obstacle. Region Vic 
lies within region VIb and has height and length 0(e5). 


and the 0(5 1 / 2 ) interior flow equations are 

Qoa = -KVPi b , 0 = V • Q 0a in region V. (2.64a, b) 

In region V, the first-correction to the boundary condition on the channel walls is given 

by 

W 0a = 0 on x = 0,1 for X + < 0, (2.65) 

while the far-field conditions are obtained by matching with region IV, to yield 

Pib —> 0, Q 0a —> 0 as X + —> —oo for 0 < z < 1. (2.66a, b) 


In region VI, the boundary conditions on the channel walls are obtained by matching 
with region VIb. We show in Appendix IB. 1 1 that. at this order, the flow problem within 
region VIb is equiva lent to the standard Prandtl problem of uniform flow past a flat 
plate ( Prandtll 190 4). The matching conditions are therefore given by 


w 0a = ft (4X+) 1/2 on 2 = 0 for X + > 0, (2.67a) 

woa = —ft (4A' + ) 1 on z = 1 for X + > 0, (2.676) 


where /3 (~ 1.721) is a constant that is determined numerically from a nonlinear ordi¬ 
nary differen t ial eq uation (discussed further in Appendix IB.II) ; the parameter /3 used 
Van Dyke] ( 19701 ) and Wilson (1971) is related to ft via ft = 2 1 / 2 /?. Finally, the inter¬ 


facial conditions onI + = 0 are 


U 0 a = u 0a , Pia = 0, w 0a = 0 on X + = 0 for 0 < z < 1. (2.68a - c) 

We do not solve the coupled system of equations (12.64[) - (12.68l) here, though we note 
that it is possible to obtain an implicit representation for wo a via a Fourier transforma¬ 
tion. Instead, we acquire the necessary information for the matching with region Via, 
which allows us to determine II(j" in the next section. We note that (12.631) can be rear¬ 
ranged to deduce the following three results: firstly, that 

V 2 woa = 0 in region VI; (2.69) 

secondly, that p\ a is the harmonic conjugate of woa', and, thirdly, that pi a +uia is indepen¬ 
dent of X + . We deduce that the leading-order far-field behaviour of Woa by considering 
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the governing equation (12.6911 together with the channel wall boundary conditions (12.67 all 
and (12. 67 foil as X + —>• oo, to obtain the expansion 

W 0a ~ —tt! as X + —> oo for 0 < 2 < 1. (2.70) 

(4X+) 1/2 

Using the far-field condition (12.701) in the governing equations (I2.63h .c'l. we further deduce 
that 

u 0 a ~/3(4X+) 1/2 , 4X+) 1/2 as X+ -> oo for 0 < * < 1. (2.71a, 6) 

We note that the far-field expansions (12.701) (12.711) have no dependence on the per¬ 
meability K. As these far-held expansions are used to form matching conditions, which 
transfer information into the transition region Via (which is where we calculate the pres¬ 
sure drop IIq ), we note that IIJ will be independent of K. We carry out this matching 
in Appendix IB. 21 

The components of stress on the interface at x = 1 are given by 

<t xx = - ^12 - + e6~ 1/2 p la ^j + o(e6~ 1/2 ), (2.72 a) 

<?xz = e5 1/2 ( u 0az + w 0aX +) + o(eS 1/2 ) (2.726) 

as e,5 —> 0 with e C 5 C 1, where we have not explicitly calculated p la , uo a , or w o a - 
In terms of the first order at which the stress varies in z, we find that the normal 
stress (12.72al) for outflow is a factor of 0(5 ~ 1/ ’ 2 ) larger than for inflow (12.60al) . while the 
shear stress (12.7261) for outflow is a factor of 0(5 1 ^ 2 ) smaller than for inflow (12.6061) . 

2.5.7. Transition region Via 

We now consider the region in which the flow transitions from plug to Poiseuille flow. 
This occurs when the growing Prandtl boundary layers with thickness of 0((5X + ) 1 / 2 ) 
become of 0(1), that is, when X + = 0(1/5), bringing about a balance in the exterior flow 
equations between fluid momentum and viscosity across the full channel width. Recall 
that we have called this region Via, as illustrated in figure |TT] 

The relevant scalings are X + = 5~ 1 X and Wo = 5 wq, and we introduce the new 
variables u o = u$ and p\ = p\ to signify that we are working in region Via. Then, in 
region Via (where 0 < X < oo and 0 < 2 < 1) the exterior flow equations given in 
figure [5] become 

Mo u 0 x + w 0 U0z = -dPix + $ 2 %xx + Mo Z z, (2.73a) 

6 2 (u 0 u 0 x + w 0 w 0z ) = -6p lz + 5 4 w 0 ^ + 6 2 w 0zz , (2.736) 

m 0 _y + wo z = 0. (2.73c) 

We pose the asymptotic expansions 

mo = woo + 0(5 1/2 ), w 0 = w 0 i + 0(5 1/2 ), pi = 5 _1 pio + 0(5 _1/2 ), (2.74a - c) 

as 5 —> 0, and substitute them into the governing equations (12.731) to yield the following 
leading-order governing equations in region Via: 

MooMqqJ^ Woi'Uqo z P\0X MOO zzi 0 PlOzi M qox woiz 0- (2.75a c) 

The boundary conditions on the channel walls become 


Moo = 0 on z = 0,1 for X > 0. 


(2.76) 
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The leading-order matching conditions are derived in Appendix IB.2I They are ob¬ 
tained by forming com posite expa nsions between regions VI and VIb, then using Van 
Dyke’s matching rule ( Van Dvkel 1975 ). We use the multiplicative composite expan¬ 
sions ( Van Dvkgll975 ) to ensure that the matching condition for the velocity satisfies 
the no-slip condition on the channel walls. This allows a greater accuracy when we solve 
numerically the resulting system. The matching conditions are given by 


woo ~ (l + 2/3X 1/2 ) /' (t7 (0) ) /' (v( i)) , 

-or ~ ^1/2 (mf (^(0)) - / (’7(0))) (/ (7(1)) - 7(i)/' (7(1))) , 

Pio ~ -2/3X 1 / 2 , 


(2.77 a) 
(2.776) 
(2.77c) 


as X l 0, where / is the standard Blasius similarity solution (defined in (IB 71) 1. p( 0 ) = 
z/X 1 ! 2 and pm = (1 — z)IX x l 2 . The matching conditions (12.771) are consistent with a 
small X expansion (for fixed *) of the governing equations (12.751) . The matching condi¬ 
tions between regions Via and VII are given by 

woo ~ 6z(l—z), w 01 —> 0, pio ~ —12A’’+IIo" as X —> oo for 0 < z < 1, (2.78a — c) 


where IlJ is determined from the solution in region Via. The full system for the flow in 
regi on Via is given by (12.751) (12.781) . 

In Bodoia fe Osterkl ( 1961 ), this same transition region was considered, though was not 
derived formally through asymptotic methods. The authors used the boundary conditions 


woo = 1, woi = 0, pw = 0 on X = 0 for 0 < z < 1, (2.79a — c) 

instead of the matching conditions (12.771) . We note that, unlike (12.7761) . the condi¬ 
tion (12. 79b ) is not consistent with the small X approximation to the governing equa¬ 
tions (12.751) . from which it can be deduced that woi is singular at t h e corners at X = 0 
z = 0,1. With the boundary conditions (12.791) . it was found in iBodoia fe Osterlel (1961) 
that the pressure jump II^ « —0.338. 

We now calculate IIq~ numerically using the formal matching conditions (12.771) . We 
start the numerical simulation at X = 10~ 4 to avoid the inverse square-ro ot sin gul arities 
in woi and in p 10 y> and use the finite-difference scheme outlined in B odoia fe Osterle 
( 196ll) . In figure fl2l we show the numerical solutions for pio + 12X as a function of X, 
and we emphasise that p io is independent of z. From the matching condition (I2.78fc h 
we determine that the value of the leading-order outflow pressure jump between outer 
regions IV and VII is given b y Il f t = —0.3 31 to 3 significant figures. We note that 
the value of II + found bv iBodoia fc Osterlel (1961) using their inconsistent boundary 
conditions (12.791) only differs by about 2% from the value of II + obtained using formal 
matching conditions. 

Combining the original asymptotic expansion (12.261) with the asymptotic expansion in 
region Via (12.74b ). we deduce that the pressure jump between outer outflow regions IV 
and VII is 0(e/S). Specifically, we have determined that 


n+ ~? 

o 


as e, S —> 0 with e -C <5 <?C 1, 


(2.80) 


where IlJ 


-0.331. 
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Figure 12. Calculating the outflow pressure drop along the channel by plotting pio + 12A. 

See text for further details. 

2.6. Summary 

We have thus far determined the boundary-layer structure for two-dimensional flow 
through a channel with stationary walls containing a porous obstacle. We formulated the 
full problem for an 0(1) Reynolds number, then considered the asymptotic sub-limits in 
which e <C Re -C 1, and 1 -C Re 1/e. In the first sub-limit, the inflow and outflow 
boundary-layer structure was found to be the same, whereas in the latter sub-limit, a 
rich boundary-layer structure was unveiled, which differed for inflow and outflow. We 
determined how flow transitions from Poiseuille to plug for inflow and from plug back to 
Poiseuille for outflow, which allowed us to determine the leading-order asymptotic values 
of the pressure jumps II and II + across regions II/III and regions V/VI, respectively. 
For inflow, where the flow problems in regions II and III decouple, we calculated the stress 
acting on the interface up to the order at which there is a variation in z. For outflow, 
where the relevant flow problems in regions V and VI are coupled, we determined the 
scalings of the stress acting on the interface up to the order at which there is a variation 
in z. 

In the next section, we generalize the problem and consider unsteady three-dimensional 
flow in a Hele-Shaw cell past a tight-fitting, highly permeable cylindrical porous obstacle, 
with a smooth boundary. As the fluid can now travel around the porous obstacle, the 
high-Reynolds-number limit has an entrainment effect for outflow. We show that this 
entrainment can lead to a net force acting on the obstacle, even when the unidirectional 
far-field forcing is periodic with a zero mean. 


3. Unsteady three-dimensional flow 

A schematic of the three-dimensional problem is illustrated in figure [2 ]d, and a two- 
dimensional plan view is illustrated in figure 1131 We initially work in a Cartesian co¬ 
ordinate system (a;, y, z) with origin at the obstacle centre of mass projected onto the 
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Figure 13. Plan view of streamlines through and around a general porous obstacle model with a 
uniform far-held how. The 2 -axis points out of the page. The Hele-Shaw cell has a depth of length 
h in the 2 -direction (not shown). We show the unit vectors e x and e y in the x- and y-directions 
respectively (for the Cartesian coordinate system), and the curvilinear (n, s)-coordinate system 
for a given point on the boundary of the obstacle. 


planar plate of the Hele-Shaw cell in the plane 2 = 0, the other plate lying in the plane 
z = h. The obstacle axis is normal to the two parallel flat sides of the Hele-Shaw cell, 
and is parallel to the 2-axis. The typical cross-sectional extent of the obstacle is l and 
the Hele-Shaw cell height is h , and we shall also assume that £ = /i/l<l. 

The exterior and interior fluid velocities are denoted by u s = u s e x + v s e y + w s e z and 
Q s = U s e x + V s e y + W s e z , respectively. The exterior and interior fluid pressures are de¬ 
noted by p s and P s , respectively. We use the superscript s to denote variables in the three- 
dimensional problem, and we later reduce the three-dimensional problem by relating these 
variables to their two-dimensional counterparts in the three-dimensional analogue of the 
inner regions considered in 1 )2.3112.51 The typical magnitude of the unidirectional far- 
held how is Uoo- We nondimensionalize by scaling the variables as follows: (x,y) ~ l, 
z ~ el, t ~ l/Uoo (corresponding to an 0(1) Strouhal number), ( p s ,P s ) ~ //0’ oo /(e/i), 
(u s , V s , U s , V s ) ~ Uoo, and (’U> S ,W S ) ~ eUoo ■ With these scalings, the dimensionless 
Navier-Stokes equations, which hold inside the Hele-Shaw cell but outside the porous 
obstacle, are given by 


cRe (• u s t + u s • Vu s ) = -p% + e 2 Viu s + u s zz , (3.1a) 

eWe « + u s • V« s ) = -p s y + e 2 \7 2 ± v s + v s zz , (3.16) 

e 3 Re (w% + u s • Vw s ) = —p s z + e 4 V^u> s + e 2 w s zz , (3.1c) 

0 = V-m s , (3.Id) 


where Re = phUoo/P- is th e global Reynolds number, V is the three-dimensional gra¬ 
dient operator, and Vj_ = d xx + dyy is the two-dimensional Laplacian operator. The 
dimensionless Darcy equations, which hold inside the porous obstacle, are given by 

U s = —KPf, V s = —KP y , e 2 W s = —KP Z , 0 = V • Q s , (3.2a - d) 


where K = k/h 2 -C 1, as before. 
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3.1. Boundary conditions 

The boundary conditions on the plates of the Hele-Shaw cell at z = 0, 1 are those of no 
slip and no flux for the exterior flow and of no flux for the interior flow, so that 

u s — 0 onz = 0,l outside the obstacle, (3.3 a) 

W s = 0 on z = 0,1 inside the obstacle. (3.36) 

On the interface, we impose three-dimensional equivalents of the two-dimensional in¬ 
terfacial conditions given by m- These are continuity of normal flux, continuity of 
pressure, and a no-tangential-slip condition, as follows 

u s • n = Q s • n on the interface, (3.4a) 

pS _ ps Qn j n t er f aC e 5 (3.46) 

u s — ( u s • n) n = 0 on the interface. (3.4c) 

where n is the unit normal pointing out of the porous obstacle. 

The far-field condition is 

p s ~ —12 xG(t) as x 2 + y 2 —> oo, (3.5) 

where G(t) is a dimensionless function of time. We take |G(t)| and |G'(t)| to be of 0(1), 
so that the far-field forcing does not have a large velocity or acceleration. 


3.2. Asymptotic structure 

The dimensionless problem ( 13 . 11 ) ( 13 . 51 ) is characterised by two lengthscales: the typical 
cross-sectional extent of the obstacle, 1, and the Hele-Shaw cell height, e 1. We again 
seek a solution using the method of matched asymptotic expansions in terms of the small 
parameter e. We note that, in contrast to (J2] not all of the fluid has to pass through the 
porous obstacle as it is now able to flow around the porous obstacle. In this section, 
we show that this difference results in a Darcy velocity of O(K) with an 0(1) pressure 
drop, in comparison to the two-dimensional case where the Darcy velocity was of 0(1) 
with an 0(1/K) pressure drop. As the continuity of flux condition (I3.4al) ensures that 
u s ■ n = O(K) near the interface, we are able to deduce that K = 0(e) corresponds 
to an obstacle which is impermeable at leading order in e. As discussed previously, we 
wish to investigate the case where the obstacle is as permeable as possible, and we 
therefore consider the limit where e -C K <C 1. The local Reynolds number (defined 
later) characterises the flow close to the interface. We initially consider the distinguished 
limit whereby the local Reynolds number is of 0(1), then consider the asymptotic sub¬ 
limits of a large and small local Reynolds number in more detail, as before. 

The outer problems are characterised by an 0(1) lengthscale in the (x, j/)-plane (with 
the Hele-Shaw cell height, e, being the lengthscale in the z-direction), where we have at 
leading order Poiscuillc and plug flow in the exterior and interior regions, respectively. We 
assume that the obstacle is a cylinder whose cross-section has a curvature of 0(1), so that 
the surface near a point on the obstacle boundary is well-approximated by the tangent 
plane at that point. The inner regions are then characterised by an 0(e) extension of each 
local normal plane to the interface and occur on either side of the interface, so that the 
inner problems are similar to those considered in <51 and the boundary layer structure is 
its three-dimensional equivalent. 
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3.3. Outer regions 

In the outer regions, we take asymptotic expansions in the limit as e —> 0, of the form 

/ = fo + e/i + 0(e 2 ), (3.6) 

for / € {u s , V s , w s , U s 7 V s , W s ,p s , P s }. At leading-order we find that, exterior to the 
obstacle, 

u s 0 = 9 i(z)p s 0x , = gx{z)p s 0y , Wq = 0, V • wg = 0, {3.7a-d) 

where gi{z) = z(z — l)/2, while interior to the obstacle, 

U$ = -KP° X , Vq = —KPq v , Wq = 0, V-Q s 0 = 0, (3.8 a-d) 

so that both the interior and exterior pressure satisfy Laplace’s equation in two di¬ 

mensions. We define the two-dimensional cross-section of the obstacle as D and the 
one-dimensional boundary of this as dD. Thus, the governing equations are 

V \p s 0 = 0 outside D, =0 in T. (3.9a, b) 

The far-held condition is 

Pq ~ —12xG(t) as x 2 + y 2 —> oo. (3.10) 

In 1 13.41 we find that the leading-order coupling conditions on the interface are 

r)n s d P s 

p s o = p o , ^ = 12 K-^j- on dD, (3.11 o, 6) 

where (13.11b ! is equivalent to the continuity of total flux condition given in (12.111 i'l. 
The leading-order problem is defined by (13.91) (13.111) . and can be solved for a given 
obstacle boundary dD. For future reference, we deduce from (13.81) that the segments 
of the boundary where we have inflow/outflow at leading order occur when dPg/dn 
evaluated on the boundary is positive/negative. 

For the O(e) correction terms we find that, outside D , 


u s i = (pi9i(z) + Re (p s ot92{z) + \S7p s 0 \ 2 g 3 {z)^ , 

(3.12 a) 

v i = §fj (pi9i(z) + Re(p s ot92{z) + |Vpg| 2 g 3 (z)')'j , 

(3.126) 

yjf = -Vi (pt J* 9 l {t) df + lte |Vpg| 2 j * g3 (0 d^ , 

(3.12c) 

V • ul = 0, 

(3.12 d) 

where 


92 {z)~ (z 2 z + i), 

(3.12e) 

!»(*) - Mo’ 4 -’ 3 + -’ 2 + = + 1 ). 

(3.12/) 

while in the obstacle, 


Ui=-KPf x , Vf = -KPf y , W/ = 0, V-Q\= 0 in D. 

(3.13a — d) 

so that the pressures satisfy 


Vi (rf + IVpSI 2 ) = 0 outside D, V 2 ± Pf = 0 in D. 

(3.14a, 6) 
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The time derivative of V 2 pg does not appear in (13.14k l due to the governing equation 
(13.9k l. The far-held condition is 

p{ —> 0 as x 2 + y 2 —> oo. (3.15) 


The first-order coupling conditions on the interface are determined in T3.4I and Ap¬ 
pendix [Cj and are given by 


p{ - PI = n s ~H(-u 0 • n) A n S+ H{u 0 ■ n), 


d ( . Re 

d^{ p '-To 



- 12 K 


dP( 

dn 


= A s H(—UQ-n) + A s+ H(uQ-n), 


(3.16a) 


(3.166) 


for Uq • n ^ 0, where a superscript ending in —/+ refers to an inflow/outflow quantity, 
and H(x) is the Heaviside step function. These inflow/outflow segments of dD are deter¬ 
mined from (13.711 (13.1111 as previously discussed. The presence of the Heaviside functions 
in (13.161) introduces the possibility of pressure singularities appearing at points on the 
boundary at which Uq • n = 0. We discuss this in more detail once we determine the 
asymptotic approximations of the coefficients pre-multiplying these functions. 

The coupling condition (13.1661) represents conservation of mass across the interface at 
this order and, in contrast to the two-dimensional versions (12.19k ~l. (12.20h ~). has addi¬ 
tional terms which correspond to an entrainment effect at this order. That is, the jump 
in the 0(e) average normal velocity across the interface (between outer regions) has 
a contribution from the variation of the leading-order tangential velocity in the inner 
regions. 


3.4. Inner regions 

As the inner regions are close to the obstacle boundary, it will be convenient to work in 
general curvilinear coordinates (n, s, z ) such that n = 0 on the boundary, with n > 0 
corresponding to the exterior of D , and s being the arc length along the boundary 
measured anticlockwise, as illustrated in figure [13] 

The relevant inner scalings are n = eN and (w s , W s ) = ( w s , W s )/e, and we define the 
components of fluid velocity in curvilinear coordinates by u s = u s n+v s t+w s e z and Q s = 
U s n + V s t + W s e z in the exterior and interior regions, res pectivel y . Under these scaling s, 
the exterior-flow equations (13.11) become (see, for example, Schlichting &; Gerstenl ( 2000 )1 


Re ( u Uj y A w u 2 ) A Oi^eRe) — e p^ A 'Unn (3.17a) 

Re ( u Vjy A w u 2 ) A 0(6i?e) = — p s A 4" v s zz A O(e), (3.176) 

Re ( u s w s N A w s w s z ) A O(eRe) = —e~ 1 p s z A w s nn + + O(e), (3.17c) 

0 = u s N A w! A e ( k(s)u s A 0*) A 0(e 2 ), (3.17d) 

where k(s) denotes the curvature of the obstacle boundary (positive if the centre of the 
osculating circle lies in the region in which n < 0 for a given s), and we assume that 
|k| = 0(1). Similarly, the interior-flow equations (13.21) become 

eU s = -KP %, (3.18a) 

V s = -KP S S A O(e), (3.186) 

eW s = ~KP S Z , (3.18c) 

0 = U S N A W s z A e ( k(s)U s A V s s ) A 0(e 2 ). (3.18d) 
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The boundary conditions on the Hele-Shaw cell walls (13.311 become 

u s = 0 on ^ = 0,1 for N > 0, (3.19a) 

W s = 0 on z = 0,1 for N <0, (3.196) 

The interfacial conditions (13.41) on N = 0 become 


u s = U s , p s = P s , v s = 0, W s = 0. 


(3.20a — d) 


We form inner expansions in powers of e as follows: 

(u s ,Q s ,p s ,P s ) = (BS.eS.tf.P'o) + e(ul,QlPl,Pl) + 0(e 2 ) as e -> 0; 


(3.21) 


substituting these expansions into the inner exterior-flow equations (I3.17h .cl and the 
interior-flow equations (I3.18h .cl. then equating powers of e, yields the leading-order pres¬ 
sure equations 


0 — ~P0Ni 


o = -i4, o = -p s on , o = -Pj 


0 z- 


(3.22 a - d) 


This leading-order system is the same as in the two-dimensional case, and we proceed 
in exactly the same way. Matching with the outer pressures, and using the leading-order 
version of the continuity of pressure condition (13.20b ). we deduce that the leading-order 
pressures are unchanged through the inner regions and we are justified in writing the cou¬ 
pling condition (13.11h l. In a similar manner, integrating the leading-order version of the 
continuity equations (I3.17J) and (I3.18dl) over the cell height, using the cell-wall bound¬ 
ary conditions (13.191) . and then matching with the outer regions, justifies the coupling 
condition (I3.11b l. The entrainment effect does not occur at this order. 

The first-order flow equations in the (N, z)-plane decouple from the tangential veloc¬ 
ities v s and V s , and the tangential coordinate s can be treated as a parameter in the 
problem. We can further reduce these problems to the two-dimensional cases considered 
previously using the scalings 


(uo(N, s,z,t),U s (N, s,z,t)) = ~v(s,t) (u 0 (N,z),U 0 {N,z )), (3.23 a) 

(w s 0 (N, s, z, t), Wq(N, s, z, t)) = \v(s, t) | (w 0 (N, z),W 0 (N , z)) , (3.236) 

(. P s i(N,s,z,t),Pl(N,s,z,t )) = (Pi(0,s,t),Pi(0,s,6)) + Ks,t)| (pi(iV, 2 :),Pi(IV,z)) , 

(3.23c) 


where the terms on the right-hand side without a superscript s are variables from the 
two-dimensional problem, and v(s,t) = - J^Mo • ndz = pg„(0, s, t)/12 = KPg n (0,s,t) 
is positive for inflow and negative for outflow. With N = —X~ for inflow and N = X + 
for outflow we obtain the systems shown in figures |4] and [5l but with Re replaced by 
the appropriate local Reynolds number, namely Re(s,t) = l^jPe. The three-dimensional 
pressure jump functions in (13.16al) are related to their two-dimensional equivalents via 
the expression 

(n s -(M),n s+ (M)) = KM)I (n-,n+). ( 3 . 24 ) 

We have, therefore, already calculated the asymptotic behaviour of II s- and II S+ for the 
pressure coupling condition (13. 16 all in 1)2.5.11 and 1 )2.5.51 We see that n s_ ,II s+ — > 0 as 
v —> 0 and thus (I3.16al) does not impose a pressure singularity at points where v = 0. 
However, the possibility still remains that (13.1661) could still impose a pressure singularity 
when v = 0, in which case an additional inner region would be present. We proceed by 
assuming that \u\ = 0(1), and discuss any pressure singularities when they occur. We 
now calculate the behaviour of the tangential velocities Vq and Vq within the inner regions 
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(which are new to the three-dimensional problem). This will allow us to determine local 
properties near the interface, for example, the stress, and to couple the outer regions by 
determining the functions A s_ and A s+ in the final coupling condition (13.1661) . 

Within the obstacle, the tangential velocity is given at leading order by 

V s o = -KP^(0,s,t). (3.25) 

The scaling 

Vo{N, s, z, t) = -04,(0, s, t)/!2)v 0 (N , z; s, t), (3.26) 

where we reiterate that s and t are now parameters in the inner three-dimensional prob¬ 
lem, results in the exterior-flow equation 

Re(s , t) (u 0 v 0 x + wovoz) = 12 + v 0 xx + v 0zz . (3.27) 

Here, we introduce X such that X = X~ < 0 for inflow, and X = X + > 0 for outflow, 
thus remaining consistent with the two-dimensional analysis. We emphasise that, as with 
the two-dimensional case, we expect the behaviour of the tangential velocity to vary 
depending on whether we have inflow or outflow. The cell-wall boundary conditions for 
Vo are given by 

v 0 = 0 on z = 0,1 for X~ < 0 and A+ > 0, (3.28) 

the interfacial conditions are 


no = 0 on X = 0 for 0 < z < 1, 


(3.29) 


and the matching conditions are 

v 0 —> 6z(l — z) as A' - —> —oo and X + —> oo for 0 < z < 1. (3.30) 


In Appendix [C] we integrate the O(e) continuity conditions (I3.17cij) . (I3.18dll over the 
channel height to obtain the expression 


d_ 

dn 


(pi - f (* -1 ivrfl 2 ) 

/>00 n 1 

= 12 / (v s 0s (N,s,z,t)~ 

Jo Jo 



dP[ 

dn 


u s Qs (0,s,z,t) ■ s) dz dN , 


(3.31) 


encompassing the entrainment effect which occurs at this order. Using the outer solu¬ 
tions m and the tangential velocity scaling (13.261) in (13.311) , we deduce that the inflow 
and outflow flux jumps A s_ and A s+ in (13.1661) are given by 


o 1 

A s ~ = - J J-^-(p s Os (0,s,t)(y o (X-,z;s,t)-l)) dzdX~, (3.32 a) 

— oo 0 
oo 1 

A S+ = -f J j- s {p s os(0,s,t)(v 0 (X + ,z-,s,t)-l)) dzdX+, (3.326) 

0 0 


where vq satisfies the system of equations (13.271) (13.301) . We proceed by determining the 
leading-order asymptotic behaviour of A s ~ and A s+ in the sub-limits in which Re <C 1 
and 1 < 1/lf < & C 1/e. 
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3.5. Small local Reynolds number: Re -C 1 

For the small local Reynolds number case, Re -C 1, the inertial terms in the governing 
equation (13.271) are not present at leading-order in Re. Therefore, no has no dependence on 
s or t (which come from Re(s,t)), nor on whether we have inflow or outflow. The local 
Reynolds number uses the local velocity component normal to the interface averaged 
over the channel height. Thus, the problem we present in this section governs the case 
when the normal local velocity is small, as well as the case when the global Reynolds 
number is small. The former case can occur, for example, in the transition between 
inflow and outflow for a large global Reynolds number, or when the obstacle is close 
to impermeable. For the l atter reason , this problem also appears in the impermeable 
obstacle case considered in iThompsonl ( 19681) . 

We consider outflow without loss of generality, the 0(1) terms in (13.271) yielding the 
governing equation 


-12 = v 0 xx + v 0zz , with X+ = X > 0, 


(3.33) 


It follows from (13.331) and the outflow versions of the boundary conditions given in (13.281) 
(13.301) that 


vo 


= 62(1 — z) — 24 a k 3 exp(—afcA') sin 0 ^ 2 , 


(3.34) 


k =0 


where ak = (2 k + l) 7 r for non-negative integers k. Using this flow result, we can deduce 
that the shear stress in the tangential direction is 


ev s ox (0,z-,s,t) = e -2p s 0s (0, s, t) ^ a k 2 sina fe 2 


(3.35) 


k -0 


as e, Re —)> 0. 

Substituting the tangential flow solution (13.34k ) into the outer flux jump conditions 
(13.321) . we obtain at leading order an analytic expression for the coupling condition 
(13.16&I) . namely 


dpi 

dn 


- 12 K 


dP( 93C(5) 


dn 


27r 5 


-Pos 


on dD , 


(3.36) 


where is the Riemann zeta function. As inflow and outflow are reversible in Stokes 
equations, the Heaviside functions in (13.16&I) disappear and we can deduce that, in the 
limit in which e -C Re <C 1, no pressure singularities exist up to 0(e). 


3.6. Large local Reynolds number: l<Cl/7f<CRe<ICl/e 

Now we consider the local asymptotic sub-limit 1 -C 1 /K <C Re -C 1/e and, as in 112.51 
we find it useful to work with the inverse local Reynolds number, S s = Re~ x -C 1, for 
ease of notation. The asymptotic structure here is the same as in 12.51 our goal is to 
determine the behaviour of the tangential velocity vo from the system (13.271) (13.301) , and 
use this to determine the asymptotic behaviour of the flux jump functions from (13.321) . 

We scale the first correction terms in the outer regions with the global Reynolds number 
as follows 

(ul,Q s 1} pl,P?) = Re(ul 0 ,Ql 0 ,pi 0 ,Pl 0 ) + (ul 1 ,Ql 1 ,pl 1 ,Pt 1 ) + 0(Re~% (3.37) 

but the inner regions will be asymptotic expansions in the local Reynolds number. We 
split the analysis into inflow and outflow, starting with the former. 
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3.6.1. Inflow 

The inflow behaviour is similar to the flow considered in 42.5.11 for which the asymp¬ 
totic structure is given in figure [7] We show that the leading-order tangential velocity is 
independent of X~ through region II and region III, but the next-order flow depends on 
X~ in region II. In this section we also show that the no-tangential-slip condition (13.291) 
is satisfied in region Ha, and we calculate the asymptotic behaviour of A s ~ as S s —»• 0. 

In region II, the tangential flow is unchanged at leading order and satisfies the asymp¬ 
totic expansion 

vo(X~,z; s,t) = 6z(l — z) + 5 s voi + 0(5^ 3 ) as S s — > 0. (3.38) 

Substituting (13.381) into the governing equation (13.271) and equating terms of 0(1), gives 

z{ 1 - z)v oix- + (1 - 2z)wqi = 0, (3.39) 

with the matching condition 

v 0 i —y 0 as X~ — > — oo. (3.40) 

Equation (13.391) . with far-held condition (13.401) . can be integrated numerically with re¬ 
spect to X~ to determine uoi- The resulting correction to the tangential velocity induces 
a slip on the cell walls, which is remedied by a boundary layer in region lib (as illustrated 
in figure [7]), the details of which are omitted from this paper for the sake of brevity. 

The coordinate scaling in region Ha is X~ = 5 S X, and we use hats to denote all 
variables in this region. We use the how variables calculated in 1 12.5.31 and pose an 
asymptotic expansion for the tangential velocity of the form 

v 0 (X,z;s,t) =v 0 o (X,z-,s,t) + 0(S s ), (3.41) 


where the leading-order term satishes the governing equation 

6 z(l - z)v 0Q x = v 00 xx for X < 0, 0 < z < 1, (3.42) 

with no-tangential-slip condition (13.291) yielding 

Voo = 0 on X = 0, 0 < 2 < 1, (3.43) 


and a matching condition with region II given by 

voo(X, z\ s, t) ~ 6z(l — z) as X — > — oo, for 0 < z < 1. (3.44) 

The system (13.421) (13.441) is solved by 

vqo = 62(1 — z) ^1 — exp( 6 z(l — z)X )^ . (3.45) 


The leading-order contribution to A s ~ in (13.32 al) is of 0(S s ), and comes from both the 
0(S s ) term in region II and the 0(1) term in region Ha, as region Ha is 0(S s ) smaller than 
region II. Therefore, substituting the asymptotic expansion (13.381) into (13.32a[) . and using 
the governing equation (13.391) to relate the velocity V 01 to w 01 (which has already been 
calculated numerically in 42.5.21) . along with the leading-order solution for the tangential 
velocity in region Ha (13.451) . the leading-order contribution to (13.32 al) becomes 


A s ~ - 


I r, 

Re ds 


l , J 

?4(o,s,t) 


V 




1 - 


V 


—00 0 


2z-l 
2(1 ~Z) 


X~ 


woi(£, z\ s, t)d£ dzdX 


(3.46) 
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as e,l/Re —> 0 with e < 1/& C it C 1. We note that the integral is finite because 
Wqi — > 0 as z — > 0 ,1. 

In contrast to the two-dimensional case, the pressure now varies in the s-direction. 
Therefore, the component of the interfacial stress in the normal direction is different 
from the two-dimensional case (12.60 all , and in three-dimensions is given by 

Onn = - (P o s (0, s,t) + 6 (ReP* 0 (0, s, t) + p s n (0, s, z, t)) ) + 0(e/Re 1/3 ) (3.47) 

as e, 1/Re —> 0 with e <C 1/Pe CUCl. We note that Pg and Pf 0 have no ^-dependence. 
Therefore, as in two dimensions, the first order at which the normal stress varies in z 
is at 0(e). Finally, we determine the leading-order term in the component of the shear 
stress on the obstacle boundary in the s-direction for inflow, a ns . From the tangential 
flow solution in region Ha given by (13.4511 . and recalling the scaling (13. 26ft . we deduce 
that 

cr ns ~ 2 ; s , t) = -eRe (3 v(s, t)p s 0s ( 0, s, t) (z( 1 - z)f^j , (3.48) 

as e, 1/Re —t 0 with e <C 1/Re CltCl. 


3.6.2. Outflow 

Outflow occurs in a similar manner to > 12.5.51 for which the asymptotic structure is 
illustrated in figure fill The tangential velocity within region V is plug flow (as given 
by (13.251) 1 and the tangential flow in region VI is of 0(5 S ). The tangential flow then 
transitions to Poiseuille flow in region Via, where it becomes of 0(1). In this section 
we solve for the leading-order tangential flow for outflow and calculate the asymptotic 
behaviour of A s+ as 5 S —> 0. 

In region VI, we use the flow variables calculated in > 12.5.61 as well as the asymptotic 
expansion 

v 0 (X+, z-, s , t) = 5 S (12X+) + 0{5 3/2 ) as 5 S -> 0, (3.49) 

which satisfies the no-tangential-slip condition (13.291) . The no-slip condition on the cell 
walls is satisfied in region VIb (as illustrated in figure flTT) via a similarity solution, as 
described in Appendix IB. II As with the three-dimensional inflow case considered in the 
previous section, the component of the interfacial stress in the normal direction is now 

cJnn = - (Pq (0, s, t) + elfeP? 0 (0,s,t)) + 0(e/^/ 2 ), (3.50) 


as e,l/Re —> 0 with e <C 1/Re <C 1, and where PJ and Pf 0 have no ^-dependence. 
Therefore, as in two dimensions, the first order at which the normal stress varies in 2 
is at 0(e/dy 2 ). From (13.491) and the scaling (j3.26|) . we deduce that the stress on the 
interface in the tangential direction is 


Pq s ( 0.M) 

E 12 


^ 0 A '+(0 ,z\s,t) = 


'(s, t)R. 


= (PcUO’M)): 


(3.51) 


as 6,1/Re —> 0 with e <C 1/Pe -C K -C 1. 

The transition to Poiseuille flow occurs in region Via, within which the relevant coordi¬ 
nate scaling is X + = 8f 1 X. We further use T’o = Vo to indicate we are in this intermediate 
region, and use the variables introduced in > 12. 5. 71 We substitute the asymptotic expan¬ 
sion 


^0 = ^00 + 0 ( 5 l /2 ) as S s 0, 


(3.52) 
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into the governing equation (13.2711 to obtain the leading-order governing equation 


u oov 00 x + woivooz = 12 + vqozz- (3.53) 

The cell wall boundary conditions (13.281) become 

voo = 0 on 2 = 0,1 for X > 0, (3.54) 

the matching conditions with region VI are 

Uoo —1► 0 as X l 0 for 0 < 2 < 1, (3.55) 

and the far-field behaviour is determined via the matching condition (13.301) . and is given 

by 

uoo ~ 6z(l — z) as X —> oo for 0 < z < 1. (3.56) 


The system (13.531) - (13.561) details the transition from plug to Poiseuille tangential flow. 
To determine the flux jump A s+ at leading order, the system must be solved numerically. 
We note that the system (13.531) (13.561) is independent of K, and hence only needs to be 
solved once to determine A s+ at leading order. 

Accounting for the change of variable X + = being dependent on s, we find 


d_ d Xdv d 

ds ds v ds QX ’ 


(3.57) 


and it follows that the leading-order behaviour of A s+ is given by 


r\ 

A s+ - ReA a — {vp s 0s { 0, s, t )), (3.58a) 

where 

OO 1 

Aa = /y (l-«bo(^,*)) dzdX. (3.586) 

0 0 


We calculate voo numerically using a second-order central finite-difference scheme, then 
determine A a using the trapezium rule. We find that A a ss 0.125 to three significant 
figures. 

Combining the outflow flux jump (13.581) with its inflow equivalent (13.461) . and recalling 
that v = KPfi n (0,s,t), we find that the leading-order (in S ) coupling condition (13.1661) 
in the large Reynolds number asymptotic sub-limit is given by 


d_ 

dn 


- Jo p « + So |VpS|2 


- 12 K 


dPfo 

dn 


= -KAa-^ {PonPos) H ( u 0 ‘ n ) on 3D. (3.59 a) 

The leading-order pressure jump condition (13.16al) in the large Reynolds number sub¬ 
limit is given by 


Pw - Pi, o = (^oJ 2 n+i7(«o • n) on dD. (3.596) 


We see that the coefficient pre-multiplying the Heaviside function in (13.59 al) does not 
necessarily vanish when Uq • n = 0, whereas the coefficient pre-multiplying the Heaviside 
function in (13.5961) does vanish. Thus, at these points the pressure is continuous but the 
pressure derivative along the boundary will have a log singularity. 
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3.7. Summary 

We now summarise the main results for the three-dimensional case. The general outer 
problem is governed by the two-dimensional linear equations: outside D 


Vi p s o = o, Vi (p{ + |||VpSI 2 ) = 0; (3.60a, b) 

and inside D 

ViP o s =0, VfPf = 0. (3.61a, 6) 

The coupling conditions (in outer variables) on the boundary of the obstacle cross-section 
dD are given by 


d_ 

dn 


Po ~ Po = 0, 
Pi -P? = K 


dPo_ 

dn 


II(s, t; Re, K ), 


dPo 

dn 


- 12 


Pi 


Re 

To 


* - S |Vs 


s 12 


— 12 K 


dn 

dP( 

dn 


= 0, 


= A(s, t; Re, K). 


(3.62 a) 
(3.62 b) 

(3.62c) 

(3.62d) 


Here, d/dn denotes the outward normal derivative to the obstacle, and we define the 
functions 


U(s, t; Re, K ) = n~(s, t; Re, K)H(P$ n ) + n+(s, t; Re, K)H(-P£ n ), (3.63a) 

A(s, t; lie, K) = A S "(s, f, lie, K)H(P° n ) + A s+ (s, t; lie, K)H{-P^ n ), (3.63 b) 

where n _ and n + are outputs to the nonlinear systems presented in figures 0] and [5] 
(using Re instead of Re), A s ~ and A s+ are solutions to (13.32 aD and (13.32 61) respectively, 
and H(x) denotes the Heaviside step function. 

For e, Re —>■ 0, we found in H3.5l that = — n + and these functions can be numerically 
determined by solving a system of linear equations (shown in figure (HDa)). Additionally, 
from (13.361) . we found that 

A s ~,A s+ ~ 93 2 f P U0, S ,t) as e, Re —► 0, (3.64) 

where f is the Riemann zeta function. 

For e, l/i?e —> 0 with e<C 1/Ke<if <C 1, we found in H2. 5. II and H2.5.51 respectively, 
that 

n -~a/K, n+ ~ ReK\P^ n (0,s,t)\U^, (3.65a, 6) 

recalling that Re(s,t) = K\PQ n \Re, where a ~ —0.117 and w —0.331. In the same 
limit, we also found that A s_ = 0(1/ (ReK)) (the prefactor at this order of magnitude 
is given in (13.461) ), and that 

A s+ - -ReKA a -^- (P° n (0,s,t)p s Os (0,s,t)) (3.66) 

as e, l/i?e -» 0 with e <C l/i?e -C K -C 1, where A a w 0.125 (defined in (j3.58&|l ). 

In the limit in which e, Re 0, we found that the interfacial stress <j(ti)\qb = 
— P s qTi\o£) + O(e), where n is the unit normal vector on the interface directed towards 
the exterior fluid. The prefactors in the error term can be determined by solving a linear 
system of equations. For these stress components in the n- and z-directions, we must 
solve the problem of Stokes flow coupled to Darcy flow, and present the system to be 
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solved numerically in ^2.41 The error term due to the stress component in the s-direction 
is given in (13.351) . For the stress components in the n- and ^-directions, we showed in fig¬ 
ures [HI b) and (c), respectively, the leading-order difference between the stress predicted 
by the boundary layer analysis and the stress predicted by a naive application of the 
outer solutions. 

In the limit in which e, 1 /Re —> 0 with e <C 1 / Re <€. K <C 1, the interfacial stress is 
given by 

cr(n)\ dD ~ - (Pq(0, s, t) + eReP? o (0, s, t)) n 

- cReK (3P^(0,s,t)p s Os (0,s,t) (z( 1 - z)f) H(P° n )t, (3.67) 

where t = e z x n is the unit tangent vector on the interface in the anticlockwise direction. 
The given stress components in the normal direction in (13.671) do not vary in z , and the 
order of magnitude of the error is different for inflow and for outflow. For inflow, the 
error is of 0(e) and the correction terms in the n- and ^-directions are given in (13.471) 
and (12.6061) . respectively. In figure [lUl we show the difference between the 0(e) stress (in 
the n- and ^-directions) from the inner flow and that predicted by a naive application of 
the outer solution. For outflow, the errors in the interfacial stress are of 0(e(ReK ) 1 ^ 2 ) in 
the ro-direction, 0(e/ (ReK)) in the s-direction, and 0(e/ (ReK) 1 ^ 2 ) in the z-direction, 
as deduced in (13.501) . (13.511) . and (12.7261) . respectively. 


4. Application of results to a cylinder with a circular cross-section 
within a Hele-Shaw cell 

We now apply our results to the problem of a forced time-dependent far-held flow 
in a Hele-Shaw domain containing a porous cylinder whose cross-section is a circle with 
dimensionless radius 1. Such a set-up can arise in tissue engineering, where the interior of 
the porous obstacle is seeded with cells. Another application is in modelling the erosion 
of a porous biofilm, for example, the attempted removal of dental plaque from between 
teeth. 

In these applications, quantities of physical interest can be determined from our anal¬ 
ysis. Firstly, the net force acting on the porous obstacle determines the movement of the 
porous obstacle were it free to move. Although a periodic forcing in the l o w Rey nolds 
number regime would lead to no net movement over one oscillation Purcelll 119771) . the 
same is not true for a large Reynolds number regime, and we investigate the latter case 
here. Secondly, in tissue engineering applications, cell growth is often coupled to the shear 
stress that cells experience, which is one part of a process known as mechanotransduction. 
Thus, to help with cell placement within a porous scaffold, it is important to know how 


the internal shear stress varies within an obstacle. We use results from Whittaker et al. 


1 20091) to estimate the spatial variation of shear stress within the porous obstacle, and 
thus the spatial variation of shear stress experienced by cells placed within such an ob¬ 
stacle. Finally, biofilm erosion is ofte n taken to be propo rtional to the square root of the 
shear stress acting on the interface ( Duddu et al. 2009 ) and, to this end, we determine 
the interfacial shear stress. 

We use cylindrical polar coordinates (r, 9 , z), where r = 0 ,z = 1/2 corresponds to the 
centre of the porous obstacle, and define the components of the exterior and interior fluid 
velocities as u s = u r e r + u e eg + u z e z and Q s = U r e r + U e eg + U z e z , respectively. 

We consider the (more interesting) case where e, 1 /Re 0 with e <C 1/Re <C K <C 
1, and use the corresponding outer system as outlined in 33.71 That is, we consider 
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asymptotic expansions of the form 

u s ~ itg + eReu { 0 as e, 1/Re —> 0 with e <C 1/Pe <C K -C 1, (4.1) 

with similar expressions for Q' 9 , p s , and P s . The governing equations for the exterior 
fluid are 

Vipg = 0, Vi (pio + ^|Vpg| 2 ) =0 forr>l, 0 < 0 < 2 tt, 0 < * < 1; (4.2a, 6) 


and, for the interior flow, 

ViP 0 s = 0, ViPfo = 0 for r < 1, 0 < 9 < 2tt, 0 < 2 < 1. (4.3a, b) 


The coupling conditions on r = 1 are given by 


d_ 

dr 


P 10 


1 

10 


Po - Po = 0 , 

Pio ~ P?o = 


dp s o 

dr 


dP s 
- 12 K- 0 


* ‘ 56 ^ 


si2 


- 12 K 


dr 

dP s 

(To 

dr 


= 0, 


K dP * 


dr 


H(~PSr), 


= -K A„ 


<9(9 


(P 0 >^)i7(-P 0 s r ), 


(4.4a) 

(4.46) 


(4.4c) 

(4.4d) 


where K is the dimensionless permeability, 11))", A 0 , and At are constants defined in 13.71 
and H ( x ) is the Heaviside step function. The far-field conditions are 

Pq ~ — 12a;G(i), —>• 0 as r — > oo. (4.5) 

The leading-order problem is solved by 

Pq = —12 ^ G(t) cos0, Pg = —12(1 + A)rG(t) cos 9, (4.6a, 6) 


where A = (1 — 12A")/(1 + 12AT). The points on the obstacle boundary at which the flow 
transitions between inflow and outflow are when Pq t = 0 which, for this specific problem, 
is when 6 = 7t/2, 37t/ 2, and the region 9 £ (7r/2, 3n/2) admits inflow when G > 0, and 
outflow when G < 0. The first-order problem is solved by 

s 27 A 2 G 2 (t) ^ 

Pw = -jjjj-+ Z^ a ™(i) r cos n9, (4.7 a) 

n= 1 
oo 

Pio = ^2 b n (t)r n cos n6, (4-76) 

n =0 

where the coefficients a n , b n can be determined via a standard application of Fourier 
series, and we given b n (which are used in the following analysis) in Appendix [Pi 
We note that the even modes (apart from n = 0,2) vanish, and ( 02 , 1 + 1 , 6271 + 1 ) ~ 
(—l) n /(2n+ l) 2 (a,/3) as n — > 00 (where a and /3 are constant). Thus, at the points 
6 = 7t/2,37t/ 2 (where there is a transition between inflow and outflow) the 0{eRe) pres¬ 
sures are continuous but their derivatives with respect to 9 are singular, as expected. 
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Using (13.6711 . the total force acting on the obstacle, F, is given by 


/*2-7T 


F - 


(Pq (1, 9, t) + eReP/ 0 ( 1, 6, t)) (cos 9e x + sin 8e v ) d 9 


r 2 n 


— eReK / P(P 0 s r ) / 3 (z(l - z))~ P^l, 9, t)p s og {l, 9, t) (- sin6»e x + cos 9e y ) dzdd, 
Jo Jo 

(4.8) 


as e, 1 /Re ^ 0 with e -C 1/Pe -C K <C 1. Thus the only contribution to the total force 
acting on the obstacle from Pf 0 is from the coefficient b\(t), given by 


6 ' (t) = ~ <2 ™° + G(i)|G(t)| - (4 - 9) 
Thus, the force (14.81) is given by 

F~( 12tt(1 + A)G(t) + eRe ((1 + A) 2 G(t)\G(t)\ - nb^t)) ) e x , (4.10) 


as e, 1 /Re —> 0 with e <C 1/Pe -C K -C 1. 

If the far-field forcing is periodic, with zero mean, the force acting on the obstacle does 
not necessarily vanish. That is, if there exists T > 0 such that G(t + T) = G{t) for all t, 
and where f 0 G(t) df = 0, the net force experienced by the obstacle over one oscillation 
is given by 

j\(t) At ~ 96(1 + A) 2 (eReK) ^ ^G(t)\G(t)\ d^ e* (4.11) 


as e, 1/Pe —> 0 with e <C 1/Re CUCl. The integral on the right-hand side of (14.111) 
vanishes for single mode oscillations (for example, G(t) = cos t) but, in general, not for 
more complicated periodic functions (for example, G(t) = cost — cos 2t). Even though the 
forcing is periodic, the effect of fluid inertia is to impart a net force over one oscillation. 
This would cause the obstacle to drift were it free to move. We note that the right-hand 
side of (14.111) vanishes as K —» 0, and hence the net force would not appear at this 
order for impermeable obstacles. We can attribute the appearance of the net force to the 
entrainment effect that occurs for outflow, which we were able to capture via a formal 
boundary layer analysis. 

In tissue engineering applications, the shear stress S experienced within the porous 
obstacle is of interest, as cells are placed within the porous obstacle and their growth 
may be dependent on the shear stress. As different cells will require different levels of shear 
stress to grow optimally, it is important to under stan d h ow t he she ar str ess experienced 
by cells will vary across the porous obstacle. In Whittaker et al. 1 20091) ■ it was shown 
that the shear stress experienced by cells within the individual porous scaffold pores can 
be estimated from the Darcy velocity. In particular, the shear stress is proportional to 
the magnitude of the Darcy velocity, which can be obtained by the pressure solutions 
(14.6b ) and (14.7&I) . to obtain 


G 00 

S oc I KS7P S \ ~ 12K(1 + A)\G\ - eReK— ^(n + l)b n+1 (t)r n cos nh. (4.12) 

' ' ra =0 

We can see from (14.121) that the spatial variation of shear stress through the obstacle 
occurs at O(eReK). We show that there can be significant spatial variation of shear stress 
within the obstacle (figure ITU) . From the coefficients b n (t), given in Appendix IPl we can 
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Figure 14. The spatial variation of the Darcy velocity averaged over one temporal oscil¬ 
lation of the far-held, given by the O(eReK) term in (14.121) . This term is proportional to 
the shear stress experienced within the porous obstacle. We use far-held forcing functions of 
(a) G(t) = cos(t)/y / 7r, (b) G(t) = (cos t — cos 2f )/v27r. The functions are chosen such that 
fg W G 2 (t)dt = 1. In both hgures, we take K = 0.01, we use contour spacings of 0.05 and, as 
there are logarithmic singularities at (x, y) = (0,±1), we do not plot contours of values larger 
than 0.3. We see that the symmetry exhibited in (a) is not present in (b). This is due to the 
effect of fluid inertia. 


deduce that the spatial variation is symmetric across x = 0 when f^G\G\ d t = 0. Thus, 
all single-mode oscillations will cause an internal shear stress which is symmetric across 
x = 0 at 0(eRe), for example G(t) = (cost)/^/n (figure ITlh). but this will not be the 
case for more complicated forcing functions, such as G(t) = (cos t — cos2t)/\/27r (figure 
nab). The prefactors of these functions are chosen such that J^G 2 (t)dt = 1 . We see 
that there are apparent singularities in the internal shear stress at (x, y ) = (0, ±1), these 
are due to the logarithmic singularities that we have previously discussed. Although the 
shear stress will increase near these points, their singular nature is unphysical and the 
shear stress will be bounded; we expect that a formal realization of this could be achieved 
by considering boundary regions near these points if required. 


We can also apply the results of our boundary layer analysis to predict the erosion of 
a porous biofilm, by calculating the shear stress acting on the interface. Th e square root 


of thi s shear stress is proportional to the interfacial erosion due to the flow (|Duddu et al. 


20091) . Using the flow solution (14.61) in (13.671) . we see that the leading-order shear stress, 


r, is 


t = 216 (eReK) H(P£ r ) ((1 + A)G(t)z(l - z)) 2 sin29t, (4.13) 

where t = e z x n is the tangential vector on the interface in the direction of increasing 
9. The Heaviside step function in (14.131) means that r vanishes for 9 £ (7t/2,37t/ 2) when 
G < 0, and for 9 £ (—7t/2,7t/ 2) when G > 0. For a periodic far-field forcing with zero 
mean, the temporal average over one oscillation of the square root of the magnitude of 
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Figure 15. A contour plot of the function z(l — z) \ sin 20| 1 ^ 2 , showing the shape of erosion in the 
region (9, z) G [0,7t/2] x [0,1]. This shape is repeated in each quadrant of the circular cylinder 
interface, but the magnitude will vary, as discussed in the text. We use contour spacings of 0.05. 


shear stress is 
T 

M 1 / 2 df = Tz(l — z)\ sin20| 1 / 2 , 

r = (216 (eReK)) 1,2l -±^ 


(4.14a) 

f 0 T \G\H(G(t)) d t for 9 G (tt/2, 3tt/2), 
f 0 T \G\H{-G{t)) d t for d G (— tt/2, tt/2), 

(4.146) 


Here, (14.14al) is proportional to the initial biofilm erosion. Thus, we have deduced that 
this erosion has a shape of z( 1 — z)\ sin20| 1 / 2 , repeated in each quadrant of the interface. 
The magnitude of this erosion will be the product of T and the constant of proportionality 
linking erosion to the square root of shear stress. From the form of T, given in (14.1461) . 
we see that the magnitude of erosion is not necessarily symmetric across x = 0 and the 
asymmetry will depend on G{t). The interfacial erosion is locally maximal at (9,z) = 
(7r(l + 2n)/4, 0.5) for n G Z and, in each quadrant, decreases from these points of maxima 
with two lines of symmetry, one with constant z and the other with constant 6 (figure 

UM- 


5. Discussion 

We make extensive use of the method of matched asymptotic expansions to investigate 
the laminar flow around and through a porous obstacle in both a channel and a Hele- 
Shaw cell. In particular, we determine how the flow behaves close to an interface (in a 
region whose width is of the order of the small aspect ratio) between single-phase and 
porous flows (governed by the Navier-Stokes and Darcy equations, respectively) for both 
small and large Reynolds numbers. 

Our analysis allows us to resolve everywhere the leading-order fluid velocity, a nec¬ 
essary condition to start investigating the nutrient transport in the system, with the 
eventual goal to optimize nutrient delivery to cells within the porous obstacle. Further¬ 
more, we obtain important characteristics of the inner flow (in a sense made formal 
within the main text), such as the interfacial stress, and to determine suitable conditions 
to couple the (simpler) outer flows. The analytical approach reveals the dependence of 
the flow on the Reynolds number and the dimensionless permeability, without resorting 
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to numerically expensive parameter sweeps for an inherently nonlinear problem. Impor¬ 
tantly, we note that any variation of the interfacial stress in the direction transverse to 
the channel/cell cannot be determined solely from the outer solutions, and requires an 
in-depth boundary layer analysis. In particular, we must continue past the leading-order 
terms in the asymptotic expansions to obtain the relevant contributions to the stress. 

In the first part of this paper we consider two-dimensional steady flow through a 
channel fully blocked by a finite length porous obstacle and, in particular, examine the 
regions near the interface between single-phase and porous flow. These inner regions 
satisfy the full two-dimensional Navier-Stokes equations, and we calculate the pressure 
and flux jumps between outer regions in the asymptotic limits of small and large Reynolds 
number. 

Whilst we include results for the small-Reynolds-number limit for comparison and 
completeness, the case which is more physically relevant to tissue engineering applications 
is the large-Reynolds-number limit. In this paper, we show that the latter case also 
exhibits more interesting flow behaviour, and we unveil a rich boundary-layer structure. 
The leading-order (in terms of the large Reynolds number) transition between Poiseuille 
and plug flow occurs within the porous medium for inflow, and outside the porous medium 
for outflow. Therefore, the boundary layer structure is different for inflow and outflow, 
reflecting the irreversible nature of the large Reynolds number flow. We determined that 
the pressure jump was an order of the Reynolds number larger in magnitude for outflow 
than for inflow. There is no flux jump between outer regions, due to the restriction of 
fluid movement to a two-dimensional plane. 

In the second part of this paper we extend the two-dimensional work to consider un¬ 
steady three-dimensional flow in a Hele-Shaw cell containing a cylindrical porous obstacle 
(which fits tightly between the parallel planar plates of the cell) whose cross-sectional 
boundary is smooth. Due to the smooth boundary, the flow behaviour and boundary 
layer structure close to the interface are similar to the two-dimensional channel flow 
considered previously and we can reduce much of the three-dimensional problem to the 
two-dimensional problem considered in the first part of the paper. However, in contrast to 
the two-dimensional channel flow case, a small flux jump between the outer flows is now 
present, as the tangential flow allows a small amount of fluid to move around, instead of 
through, the porous obstacle. We summarise the main results of the three-dimensional 
problem in 33.71 including the results for the interfacial stress and the derived conditions 
required to couple the outer flows. 

In 0 we apply the general results to a three-dimensional cylinder with a circular 
cross-section in a flow with a far-held forcing. From this, we are able to make several 
physical predictions. We deduce that, due to the fluid inertia, a periodic forcing with 
zero mean could cause a drift on the cylinder were it free to move. This is relevant to 
the tissue engineering application we discussed in SI where it is important to control the 
movement of the porous obstacle, which has implications for the delivery of nutrients to 
cells within the construct. Moreover, knowledge of this long-term drift may be helpful in 
hushing blockages from pipes. We also calculate the spatial variation in internal shear 
stress through the obstacle. As cell growth can be coupled to shear stress, this has 
important implications for cell placement in tissue engineering applications. We show 
that the spatial variation in shear stress is symmetric through the obstacle for a single¬ 
mode far-held oscillatory forcing, but can break symmetry for multiple-mode far-held 
oscillatory forcings. We are also able to determine the interfacial stress and use this 
result to obtain the initial interfacial erosion of a porous biohlm, assuming that interfacial 
erosion is proportional to the square root of the modulus of shear stress with a given 
constant of proportionality. 
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In the large-Reynolds-number limit, we apply three different tangential slip conditions 
on the interface to see the effect of varying boundary condition. We find that any of the 
three tangential slip conditions we apply only produce the same leading-order coupling 
conditions for inflow. Thus, knowledge of the tangential velocity condition is important 
for inflow if local information about, say, the interfacial stress were required, and it is 
important for outflow if global information about the outer problem is required. Addi¬ 
tionally, we note that although we imposed continuity of pressure on the full problem, 
the coupling condition (13.62 al) would remain the same if we had imposed continuity of 
stress instead. This is because the viscous term in the normal stress does not contribute 
to the stress at leading order. The higher order pressure coupling condition (13.62 61) would 
remain unchanged in the large-Reynolds-number limit for the same reason, but would be 
different in the small-Reynolds-number limit. 

In this work, we have considered a porous medium within which the flow is governed by 
Darcy’s equations. We restrict ourselves to these governing equations because the tissue 
engineering application we are interested in requires a porous scaffold which maintains its 
structural integrity whilst moving within the flow. Thus, at the very least, the underlying 
solid matrix of which the por ous mediu m consists must be connected, and Brinkman’s 
equations will not apply (Au riault 20091) . However, there may be applications which are 
not as restrictive, where Brinkman’s equations may apply. In such cases, the coupling 
conditions for the outer problem would be different to those derived in this paper, and 
a new analysis would be necessary to determine these. Whilst the issue of coupling the 
plain and porous flow regions is numerically simpler (as Brinkman’s equations are able to 
be applied to both plain and porous fluid domains), the presence of the Laplacian viscous 
term means that it is more difficult to make analytic progress with the resulting partial 
differential equations. For example, the boundary layer structure may be significantly 
more complex within the porous medium. 

Although we have considered a porous obstacle whose interfacial boundary is normal 
to the channel or bioreactor walls, this may not be the case. We show in fj4] that erosion 
of an initially straight wall is not uniform in any tangent direction to the interface. More 
generally, manufacturing constraints could cause the porous insert to have non-straight 
walls. Such imperfections will only change the flow problem if the gradient of the walls is 
significant in one of the boundary layers. If this occurs, the change in problem geometry 
yields a complicated extension to the problem we have presented. However, this extension 
is vital if the full dynamics of erosion are to be considered in, for example, the shape 
evolution of a biofilm. 

The tissue engineering experiment discussed in (Jl] involves a moving porous obstacle, 
which is possible as there are small gaps between the flat sides of the porous obstacle and 
the flat sides of the bioreactor (figure [l]). Whilst our work has considered unsteady flow, 
we have restricted ourselves to a pinned obstacle in this paper, and neglected the small 
gaps. To consider the dynamic problem fully, we must investigate the role of these gaps. 
The presence of gaps would have no effect on the flow up to the asymptotic order s cons id- 
ered if these gaps were smaller than the size of the boundary layers presented ( Dalwadil 
20141 ). However, a gap height of the same order as the width of one of these boundary 


layers would significantly complicate the boundary layer problem, as a result of the signif¬ 
icant change in the geometry of one or more domains in which we solve various submodels 
arising in the boundary layer analysis. Knowledge of the work presented here allows the 
porous obstacle trajectory to be calculated once the effect of the small gaps has been 
calculated. As the obstacle movement would affect the flow, this would allow a deeper 
understanding of the dynamics of nutrient transport and, ultimately, tissue growth. 

The analytic nature of our results can significantly reduce the numerical expense of 
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solving such dynamical problems. This is because the coupling conditions and stress 
components that we have determined are all in terms of the outer variables. Hence, we 
have reduced a nonlinear three-dimensional problem to a linear two-dimensional problem. 
More generally, this work highlights how the exploitation of an underlying separation of 
scales can significantly reduce the computational complexity of a physical problem. 
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Appendix A. Additional inflow regions 

In this Appendix we resolve the additional boundary layers for inflow. We first intro¬ 
duce a boundary layer near z = 0 to resolve the issue of the 0(8) slip velocity induced 
on the channel wall in region II in 8 12.5.21 We call this boundary layer region lib. We 
use X~ = X for convenience, and note that there will be a similar boundary layer near 
z = 1. 

Before we proceed, it will be useful to analyse the fluid velocities uoi and Woi in 
region II in more detail. In the separable solutions to the governing equations (12.371) and 
(12.45al) . the fluid velocities uqi and u5oi are given by 

oo oo 

u 01 =^2C k exp(p k X~)g' k (z), w 0 i = - ^ n k C k exp(p k X~)g k (z), (Ala,6) 

k =1 1 

where the eigenfunctions gk(z) satisfy the eigenvalue problem 

9k(z) + + k£\ 9 k{z) = ° for z£ (0,1) (A2) 

together with boundary conditions gk( 0) = <?fc(l) = 0 (obtained from the no-slip chan¬ 
nel wall boundary conditions given in figure 2]). The coefficients Ck would have to be 
determined numerically using the orthogonality of the eigenfunctions. 

The relevant scalings are uo = S 1 ^ 3 u, Wq = 8 3 ^ 3 w, and z = and we introduce 

the new variable p\ = p to signify that we are working in region lib. Under these scalings, 


the bulk flow governing equations in Figure [4] become 

uu x + 5 2/3 wu z = -S 1/3 p x + uzz + 8 2/3 uxx, (A 3a) 

Suihx + S 5 / 3 wwz = —pz + Suizz + d 5 ^ 3 ibxx, (A 36) 

0 = ux + S 2 ^ 3 wz- (A 3c) 

The boundary conditions on the channel wall are 

u = 0 on Z = 0 for A < 0. (A 4) 

We use asymptotic expansions 

u ~ 6Z - 5 1/3 6Z 2 + S 2/3 u 01 + 0(6), (A 5a) 

w ~ w 0 2 + 0(5), (A56) 

p~-12X + 6 1/3 p 12 + 0(5 2/3 ) (A 5c) 
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as 8 —> 0; substituting into the bulk flow equations (1A 31) yields 


6 ( Zu mx + w 02 ) = ~Pi 2 x + u 01 zz, (A 6a) 

0 = -pi 2 z, (A 66) 

0 = Uoix + W 0 2 Z- (A 6c) 

The channel wall boundary condition (I A 41) is given by 

uqi = wq 2 = 0 on Z = 0 for X < 0, (A 7) 


and, from the general form of the flow solution in region II, EH). the matching condition 
with the horizontal velocity in region II is given by 


uoi ~^2c k exp(n k X)g k (0), as Z —>■ oo. 
fc=l 


We seek solutions of the system flA 61) (1A 81) in the form 


(A 8) 


(uoi,w 0 2 ,Pi 2 x) =^2exp(p k X) (f k (Z),-p k f k (Z),b k ), 


(A 9) 


fc =i 


which automatically satisfies (I A 6bb . (I A 6 d) . and where b k and f k {Z) are to be determined. 
Substituting (IA 91) into (I A 6al) . we obtain the differential equation 

K'{Z) + al{Zf' k {Z)-f k {Z)) = b k , 


for f k (Z), where a k = 6p k > 0, and with boundary conditions 
/*( 0) = 0, f' k ( 0)=0, f' k (oo) = C k g' k (0). 

The system (IA 101) (I A 111) is solved by 

, f n Z fo Ai ( a k v ) dv d u 


f k {Z) = C k g' k (dy- 


/ 0 °° Ai (a k s) ds 


(A 10) 
(A 11a — c ) 

(A 12) 


where Ai is the Airy function of the first kind. To determine b k , note that the general 
solution for (I A 101) is 

f k (Z ) = % + B kl Z + BMZ) + B k3 h 2 (Z), (A 13) 

where the far-field behaviour of hi and h 2 as Z —> oo is 

hi ~ Z~ 5 ^ 4 ex p ^ (afeZ) 3 ^ 2 ^ , h 2 ~ Z -5//4 exp ^ (a k Z) 3A ^ as Z —> oo. 

(.A 14 a, b) 

Applying the boundary condition (IA lib ) imposes B k i = C k g k ( 0) and B k2 = 0. The 
remaining boundary conditions (IA llk -b) yield the solution for b k and B k3 , but the 
simplest route to determine b k is to equate (IA 131) with (I A 121) as Z —> oo to deduce that 


b k = -3alC k g k (0) J J Ai (a k v) dw du, (A 15) 

0 u 

where we have used the identity / 0 °° Ai(a k v) du = (3afc) _1 . 

We are unable to impose any conditions on the interfacial boundary X = 0 within 
region lib. This issue is resolved by another boundary layer near X = 0, which we call 























45 


On the boundary layer structure near a highly permeable porous interface 

region lie. The scalings from region II to region lie are Uq = 0(6 1 / 3 ), Wo = 0(<5 4 / 3 ), 
pi = 0(6 1 / 3 ), X = 0(<5 2//3 ), and z = 0(6 1 ^ 3 ). The leading-order equations yield a 
solution of u o ~ 6ztanh(3z(— X)/8), from which w o and pi can be determined. 


Appendix B. Additional outflow regions 

B.l. Regions VIb and Vic 

In this section we solve for the flow near the outflow channel wall (to match with the 
asymptotic regions considered in two 1 H2.5.6I) and three-dimensions f ( 13.6.21) 1. in a region 
we denote region VIb. For the two-dimensional channel flow problem, this boundary layer 
is equivalen t to the classic uniform high Reynolds number flow past a flat plate considered 
bv lPrandtl ( 1904 ). We therefore consider the three-dimensional problem, which uses the 
two-dimensional results. We use A' + = X for convenience. 

The relevant scalings from region VI to region VIb are: z = 5 S z, Co = 6 s v , and 
wo = 6 S w. We also use Co = u, and pi = p to signify that we are working in region VIb. 

With these scalings, the bulk flow governing equations given in Figure [5] and (13.271) 
become 


uu x + wu s = -6 s p x + wjz + 6 s u X x , (B la) 

UVX + wVz = 12 + vgz + 5 S VXX, (B 16) 

UW X + wwz = -pz + wzs + 6 S WXX J (B lc) 

0 = ux T Wz • (Bid) 

Substituting the asymptotic expansions 


u = uo + 0(<5] /2 ), v = v 0 + 0(6l /2 ), ui = iu 0 + 0(6l /2 ), p = p w + 0(6l /2 ) 

(B2a-d) 

as 6 S —> 0, into the bulk-flow equations EH, we obtain the following leading-order 
equations 


U 0 u 0 x + w 0 U0z = -Pwx + Uosz, (B 3 a) 

uqvox + w 0 V0z = 12 + Cozz, (B 3 b) 

0 = -Pioi, (B 3c) 

0 = u O x +w 0 z- (B 3d) 

The boundary conditions on the channel wall z = 0 are given by 

u 0 = 0 on z = 0 for X > 0, (B 4) 

and the matching conditions into region VI are given by 

Co ~ 1, Co ~ 12X, pio —> 0, as z —> oo for X > 0. (B 5a — c) 

The system (IB 31) (IB 51) is solved by 

wo = f'(v), C 0 = 12 Xh(r)), iv 0 = , Pw = 0, (BQa-d) 

where p = zA' -1 / 2 is the similarity variable. The function f(rf) satisfies the following 
Blasius ODE: 

f"\r 1 ) + \f(i 1 )f"(r 1 ) = 0, /(0) = f'(0) = 0, /'(oo) = 1; (B 7) 
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the function h{r}) satisfies the following ODE 

h"{v) + - f'ivMv) = -1, h(0) = 0, h{oo) = 1. (B 8) 

The function / is well studied (see, e.g. iBovd ( 2008l ll. and its behaviour up to expo¬ 
nentially small terms as 77 —> 00 is given by 


/ ~ 77 — /3 as 77 — > 00 , 


(B 9) 


where f3 ss 1.721. This far-held behaviour, in conjunction with the solution for uiq 
in (IB 6b ). yields the matching conditions (12.671) . 

Finally, we note that the scalings from region VI to region Vic are: z ~ 6 S , X + ~ S s , 
vq ~ S 2 and pi ~ (5J 1 , thus yielding the full two-dimensional Navier-Stokes equations for 
Mo and wo, and the full version of (13.271) for vq. That is, taking Re = 1 in (13.271) . We do 
not consider region Vic further. 


B.2. Matching into region Via 

In this section we form the composite expansion between regions VI and VIb in their 
far-held as X + —> 00, and thereby obtain the matching conditions (12.771) and (13.551) with 
region Via. When discussing the abscissa in region VI, we use X + = X for convenience. 

Combining the leading-order versions of (I2.62h j and (13.381) . together with the far- 
held first-correction terms (12.701) (12.711) , we obtain the region V far-held expansions (as 
X —> 00) up to 0(Sl^ 2 ) for the how, and up to 0(S7 1 ^ 2 ) for the pressure as follows 

f t0 ~ 1 + 2/3 (<Uf) 1/2 , vo ~ 12 ( 5 S X) , vo- ^x^syl ’ Pi ~ -2p(X/S s ) 1/2 . 

(BlOa-d) 

The matching conditions are obtained as follows. We hrst form the multiplicative 
composite expansion of (IB 101) with the solution near the channel wall at z = 0 (given 
by (IB 61) ), and the equivalent solution near the channel wall at z = 1 (obtained by taking 
77 = (1 — z)(5 s X)~ 1 / 2 and reversing the sign on Wo)- Then, we write these expansions in 
terms of X = S s X, and retain the leading-order terms to yield, as X j. 0, 


Uo - 

-(1 + 2 ^) f (rj (0) ) f (rj (1) ) , 

(B 11a) 

^00 " 

- 12 Xh ( 77 ( 0 )) h ( 77 ( 1 )) , 

(B 116) 

Wo - 

- 2/3X1/ 2 ^ (0)/ ^ (0) '* f (^0))) (V(i)f (vm) /(%))). 

(Bile) 

P 10 * 

- -2 px 1/2 , 

(B lid) 

m = 

z/X 1 / 2 and rj w = (1 - z)/X 1 / 2 . 



Appendix C. Flux jump condition 

In this section we determine the flux jump in terms of the tangential velocity, as given 
in (13.311) by integrating the inner continuity equations (I3.17J) and (I3.18J) over the cell 
height. Using the asymptotic expansions (13.211) . and equating powers of e, we deduce 
that 

If u s o (N,s,z,t)dz = 0, If U s 0 {N,s,z,t)dz = 0, (Cla,b) 






































and 
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f u\ n (N, s, z, t) dz + k(s) ( Uq(N, s,z,t)dz + f Vq s (N, s, z, t) dz = 0, (C2a) 

Jo Jo Jo 


U* N (N,s,z,t)dz + k(s) / Uo(N,s,z,t)dz+ / Vq s (N, s, z, t) dz = 0, (C26) 


noting from (13.251) . the scaled solution within the porous obstacle, that Ug (N,s,z,t) = 
-KPg a {0,s,t). 

The matching conditions for Uq, Vg, ztf and U\ are given by 

z)q ~ t • itg(0, s, z, t) as N —> oo, (C3a) 

Fg - -KP£ s (0,s,t) as N ^ —oo, (C36) 

u\ ~ n • (Nug n (0, s, z, t) + uf(0, s, z, t)) as N —> oo, (C3c) 

U\ ~ n • (NQ s On (0, s, z, t) + Qi(0, s, z, t)) as N —> — oo. (C3 d) 

To formally avoid infinite terms in the limit as N — > oo when we integrate with respect 
to N in (1C 2D , it is necessary to subtract the far-field limit of (1C 21) from itself (using (1C II) 
and (1C 3D ). yielding 


jf« 


'IN 


*0n\dD 


) dz = - f {v s 0s - t-u s 0s \ aD ) dz, 
Jo 


[ (u s 1N - 

J o 


n ■ Q 


On\dD 


) dz = 0. 


(C 4a) 
(C46) 


Integrating (1C 41) with respect to N between 0 and oo (using the far-held conditions (1C 31) 1. 
we deduce that 


/»1 /»00 /»1 

/ ( n - u i\dD~ «iU=o) dz = - / / (v s 0s - t-u s 0s \ dD ) dzdN, (C5a) 

Jo Jo Jo 

f\n-Ql\ aD - Ul\ N=0 ) dz = 0 , (C5 b) 

Jo 

where outer/inner variables evaluated on dD is meant in the sense of the outer/inner 
problem. Finally, integrating the leading-order interfacial continuity of flux condition 
mm across the cell height, we can relate (1C Sal) and (1C 561) . Using the outer flow 
solutions (13.121) and (13.131) . we obtain (13.311) . 


Appendix D. Flow coefficients 

The flow coefficients b n (t), given in (14.761) . are 

27. 


b 0 (t) = -G 2 (t ) (^-(1 + A 2 ) + 36IIg AT 2 (1 + Af 

( 07 A 

— - 36K(1 + A) 2 (A a + 211+A) 

„ 6 , dG 288Jf(l + A) 2 (—l) fc / 

(1 + 12 K) b 2k+1 (t) - -(1 - A)—S ok + nm+1)2 _ 4) ( Aa 


(D la) 
(D 16) 


Sty) cwicwi, 


(Die) 

(Did) 


b 2k =0 for k ^ 0,1, 




















48 


M. P. Dalwadi, S. J. Chapman, S. L. Waters, and J. M. Oliver 


for integer values of k , and where Sij is the Kronecker delta, and all other parameters 
are defined in 0 
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